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I.  Introduction 

The  present  report  summarizes  the  work  that  has  been  carried  out  under 
the  support  of  the  AFOSR  (Grant  No.  AFOSR-85-0133)  during  the  funding 
period  Feb.  1,  1985  to  Jan.  31,  1988.  The  investigation  covered  by  the  Grant 
deals  with  the  physical  phenomena  associated  with  coupled  hydromagnetic 
waves  in  the  magnetosphere  and  the  physics  of  the  coupling  between  the  mag¬ 
netosphere  and  ionosphere.  The  research  effort  is  also  extended  to  study 
naturally  occurred  or  artificially  induced  phenomena  related  to  plasma  waves  in 
the  ionosphere  and  magnetosphere.  Our  study  focuses  primarily  on  developing 
theories  to  explain  the  observational  results,  however,  the  possibility  of  provid¬ 
ing  guidelines  to  future  experiments  is  also  considered. 

During  the  past  three  years.,  our  work  has  achieved  many  interesting 
results,  which  have  been  either  submitted  or  already  accepted  for  publication. 
Basically,  the  topics  of  investigation  are  divided  into  four  general  categories:  (a) 
cavity  modes  of  the  magnetosphere  resulting  in  the  discrete  spectrum  of  the 
resonant  ULF  waves;  (b)  a  hydromagnetic  code  for  the  numerical  study  of  the 
coupling  of  hydromagnetic  waves  in  the  dipole  model  of  the  magnetosphere; 
(c)  a  theoretical  model  developed  for  explaining  the  phenomenon  of  plasma 
line  overshoot  observed  in  the  ionospheric  HF  heating  experiments;  and  (d) 
thermal  fllamentation  instability  as  the  mechanism  for  generation  of  large-scale 
field-aligned  ionospheric  irregularities.  For  the  first  two  topics,  the  hydromag- 
nelic  wave  equations  are  analyzed  analytically  in  cylindrical  model  of  the  mag¬ 
netosphere  and  numerically  in  dipole  model  of  the  magnetosphere,  respectively. 
While  the  steady  state  eigenvalue  problem  is  studied  in  the  first  topic,  the 
second  topic  is  generalized  to  the  boundary  value  problem  considering  the  cou¬ 
pling  between  hydromagnetic  waves  in  the  realistic  geometry  of  the  magneto¬ 
sphere.  For  the  third  topic,  a  nonlinear  turbulent  theory  (resonance  broaden¬ 
ing  of  electron-wave  interaction)  is  incorporated  in  the  study  of  parameter  ins¬ 
tability  excited  by  a  powerful  HF  in  the  ionosphere.  For  the  last  topic,  the 
thermal  nonlinearity  gives  rise  to  the  mode-mode  coupling;  threshold  field  and 
the  growth  rate  of  the  instability  are  derived. 


II.  Publications 

The  following  publications  include  work  supported  by  the  present  Grant, 
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1.  S.P.  Kuo,  M.C.  Lee  and  A.  Wolf,  "A  Model  for  the  Discrete  Spectrum  of 
the  Resonant  ULF  Waves,"  EOS  Trans.,  AGU,  00(18),  347,  1985. 

2.  S.P.  Kuo,  M.C.  Lee  and  S.C.  Kuo,  "A  Theoretical  Model  for  Artificial' 
Spread  F  Echoes,”  Radio  Science,  20(3),  546-552,  1985. 

3.  M.C.  Lee  and  S.P.  Kuo,  "Simultaneous  Excitation  of  Large-Scale  Geomag¬ 
netic  Field  Fluctuations  and  Plasma  Density  Irregularities  by  Powerful 
Radio  Waves,"  Radio  Science,  20(3),  538-545,  1985. 
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4.  M.C.  Lee  and  S.P.  Kuo,  "Resonant  Electron  Diffusion  as  a  Saturation  Pro¬ 
cess  of  the  Synchrotron  Maser  Instability,"  J.  Plasma  Physics,  35,  177-184, 

1986. 

5.  S.P.  Kuo  and  M.C.  Lee,  "Parametric  Excitation  of  Whistler  Waves  by  HF 
Heater,”  Proceedings  of  Int'l.  Symposium  on  "Modification  of  the  Iono¬ 
sphere  by  Powerful  Radio  Waves,”  Suzdal,  USSR,  Sept.  1988,  pp.  109-110. 

8.  M.C.  Lee,  J.A.  Kong  and  S.P.  Kuo,  "Enhanced  Ionospheric  Modifications 
by  the  Combined  Operation  of  HF  and  VLF  Heaters,”  Proceedings  of  Int’l. 
Symposium  on  "Modification  of  the  Ionosphere  by  Powerful  Radio  Waves,” 
Suzdal,  USSR,  Sept.  1986,  pp.  146-147.  . 

7.  M.C.  Lee,  J.A.  Kong  and  S.P.  Kuo,  "On  the  Resonant  Ionospheric  Heating 
at  the  Electron  Gyrofrequency,"  Proceedings  of  Int’l.  Symposium  on 
"Modification  of  the  Ionosphere  by  Powerful  Radio  Waves,”  Suzdal,  USSR, 
Sept.  1986,  pp.  111-112. 

8.  S.P.  Kuo,  M.C.  Lee  and  A.  Wolfe,  "Dominance  of  Discrete  Magnetos- 
pheric  Cavity  Modes  on  the  Generation  of  Field-Line  Resonances,"  EOS 
Trans.,  AGU,  67(16),  348,  1986. 

9.  M.C.  Lee  and  S.P.  Kuo,  "On  the  Causes  of  Spectral  Broadening  of  Injected 
VLF  Waves,"  EOS  Trans.,  AGU,  87(16),  317,  1988. 

10.  M.H.  Whang,  S.P.  Kuo  and  M.C.  Lee,  "Thermal  Filamentation  Instability 
of  EM  Waves  in  Magneto  Plasma,"  Bull.  Amer.  Phys.  Soc.,  31(9),  1438, 

1986. 

»1.  S.P.  Kuo,  M.C.  Lee  and  F.T.  DJu*,h,  "New  Interpretation  of  Plasma-Line 
Overshoot  Phenomenon,"  USNC/URSI,  G/H2,  p.  298,  12-15,  Jan.  1987, 
Boulder,  CO. 

12.  S.P.  Kuo,  M.C.  Lee  and  F.T.  DJuth,  "Overshoot  of  HF-Enhanced  Plasma 
Lines  Due  to  Resonant  Broadening  Effects,”  IES  '87  proceedings,  249-253, 

1987. 

13.  S.P.  Kuo,  M.C.  Lee  and  A.  Wolfe,  "Spectral  Characteristics  of  Hydromag- 
netic  Waves  in  the  Magnetosphere,"  J.  Plasma  Phyics,  38,  235-243,  1987. 

14.  S.P.  Kuo,  M.C.  Lee  and  T.T.  DJuth,  "A  New  Interpretation  of  Plasma-Line 
Overshoot  Phenomena,”  Geophys.  Res.  Lett.,  14(9),  961-964,  1987. 

15.  S.P.  Kuo  and  M.C.  Lee,  "Thermal  Filamentation  Instability  Driven  by  the 
Auroral  Electrojet  Current,”  J.  Geophys.  Res.,  93(A),  265-269,  1988. 

16.  S.P.  Kuo,  M.H.  Whang  and  M.C.  Lee,  "Filamentation  Instability  of  Large 
Amplitudes  Alfven  Waves,”  J.  Geophys.  Res.,  accepted  for  publication, 

1988. 

17.  M.H.  Whang  and  S.P.  Kuo,  "The  Stability  Analysis  of  a  Finite  Difference 
Scheme  for  the  Coupled  Hydromagnetic  Wave  Equations  in  the  Dipole 
Model  of  the  Magnetosphere,"  submitted  to  J.  Computation  Phys.,  1988. 


III.  Summary  of  Work  Accomplished 

(a)  A  Model  for  the  Discrete  Spectrum  of  the  Resonant  ULF  Waves 

The  discrete  nature  of  ULF  pulsations  observed  in  the  magnetosphere  is 
theoretically  explained  by  using  a  cylindrical  model  [Ftadoski,  J.  Geomagn. 
Geoelectr.,  £2,  381,  1970).  A  second  order  wave  equation  is  derived  to  charac¬ 
terize  the  coupling  between  the  poloidal  and  toroidal  modes  of  the  hydromag- 
netic  waves.  When  the  experimentally  determined  plasma  density  radial  profile 
of  ”l/r3”  form  [Cummings  et  al.,  J.  Geophys.  Res.,  T±>  778>  i960;  Park  et  al.,  J. 
Geophys.  Res.,  8£,  3137,  1078]  is  used,  it  is  found  that  the  wave  equation  has 
two  turning  points  along  the  radial  axis.  By  contrast,  only  one  turning  point 
exists  in  the  wave  equation  fo  Kivelson  and  Southwood  [Geophys.  Res.  Lett., 
12.  49,  1985]  wherein  the  magnetosphere  is  described  by  a  box  mode)  having  a 
linearly  increasing  density  profile  with  distance.  Undamped  eigenfunctions  with 
discrete  eigenvalues  are  found  in  our  collisioniess  case.  They  are  conceptually 
analogous  to  the  case  that  has  discrete  eigenstates  when  particle  is  trapped  by  a 
potential  well  as  discussed  in  quantum  mechanics.  However,  there  is  a 
difference  between  the  two  cases.  In  the  present  case,  the  function  f(r) 
representing  the  square  of  the  wavenumber  cannot  simply  be  expressed  by  the 
difference  of  a  constant  eigenenergy  and  a  fixed  potential  function  as  in  quan¬ 
tum  mechanics.  This  is  because  the  eigenvalues  is  multiplied  by  a  spatially 
dependent  function.  Therefore,  the  potential  well  determined  by  f(r)  varies 
with  each  discrete  eigenstates. 

We  have  calculated  the  eigenperiods  of  the  most  readily  excited  modes  and 
the  most  likely  values  of  the  east-west  wavelengths.  The  wave  eigenperiods  are 
calculated  to  be  in  the  Pc3-4  range.  When  the  magnetospheric  cavity  is  con¬ 
sidered  by  a  resonant  filter,  only  discrete  parts  of  the  continuous  perturbations 
on  the  magnetopause  can  couple  to  the  localized  field  line  resonance  modes 
excited  inside  the  magnetosphere,  where  the  discrete  nature  of  the  spectrum  is 
determined  by  the  eigenmodes  of  the  cavity.  Details  of  this  work  are  enclosed 
in  the  Appendix. 

(b)  Stability  Analysis  of  a  Finite  Difference  Scheme  for  the  Coupled 

Hydromagnetic  Wave  Equations  in  the  Dipole  Model  of  the  Magnetosphere 

In  the  study  of  wave  propagation  In  Inhomogeneous  magneto  plasma,  the 
physical  system  is  generally  characterized  by  coupled  partial  differential  equa¬ 
tions.  One  specific  example  is  the  geomagnetic  micropulsation  originated  by 
the  hydromagnetic  waves  in  the  magnetosphere.  It  is  shown  that  these  waves 
in  the  poloidal  and  toroidal  modes  are  generally  governed  by  two  coupled 
second  order  partial  differential  equations.  Though  these  two  equations  can  be 
reduced  and  combined  into  a  single  second  order  ordinary  differential  equation 
for  a  simplified  cylindrical  model  of  the  magnetosphere,  there  is  no  apparent 
way  to  simplify  the  system  of  coupled  equations  for  a  more  realistic  dipole 
model  of  the  magnetosphere  in  which  a  dipole  geomagnetic  field  is  assumed. 


Moreover,  the  equation  for  the  toroidal  mode  is  shown  to  be  a  parabolic  type 
of  partial  differential  equation,  while  the  poloidal  mode  is  elliptic.  Both  equa¬ 
tions  have  variable  coefficients,  therefore,  for  practical  applications,  there  is 
considerable  interest  in  finding  or  developing  means  of  dealing  with  such  a  sys¬ 
tem  of  two  mixed-type  coupled  partial  differential  equations  with  variable 
coefficients. 

In  this  work,  a  numerical  algorithm  dealing  with  the  mixed-type  coupled 
partial  differential  equations  with  variable  coefficients  has  been  developed.  The 
stability  of  the  numerical  scheme  has  been  examined  and  solutions  for  the 
hydromagnetic  wave  equations  are  obtained  on  non-staggered  grids  by  this  algo¬ 
rithm.  The  main  contribution  of  this  algorithm  is  to  offer  an  efficient  way  to 
examine  the  stability  conditions  first  by  the  Fourier  method  and  then  check 
them  by  the  matrix  methods.  Since  these  two  methods  are  different,  one  can 
be  sure  of  the  stability  conditions.  Moreover,  this  algorithm  leads  a  way  for  a 
parametric  study  of  the  stability  conditions.  It  also  becomes  an  advantage  in 
the  sense  that  the  optimum  parameters  which  can  increase  the  accuracy  of  the 
numerical  solution  can  be  determined  together  with  the  stability  analyses. 
Using  the  developed  numerical  scheme,  the  problem  of  coupling  and  propaga¬ 
tion  of  hydromagnetic  waves  in  the  realistic  dipole  model  of  the  magnetosphere 
can  be  studied.  An  example  of  coupled  hydromagnetic  waves  in  the  dipole 
magnetosphere  has  been  considered.  Details  of  tnis  wor1'  addressed  in  the 
Appendix. 

(c)  Thermal  Filamentation  Instability  Driven  by  the  Auroral 
Electrojet  Current 

The  high-latitude  ionosphere  is  constantly  perturbed  by  the  ionospheric 
irregularities.  The  Farley-Buneman  instability  driven  by  electrojet  current  is 
known  to  be  effective  in  generating  meter  and  shorter-scale  irregularities  of 
type  I  in  the  E  region.  The  irregularities  of  type  II  are  of  slightly  larger  scale- 
sizes  and  are  generated  through  the  ExB  gradient  drift  instability.  However, 
the  growth  rate  of  the  two  stream  instability  as  well  as  of  the  gradient  drift  ins¬ 
tability  is  a  function  of  the  angle  between  the  wavevectcr  and  the  direction  of 
the  electron  drift.  The  growth  rate  becomes  zero  when  the  wavevector  is  per¬ 
pendicular  to  the  drift  direction.  Therefore,  it  is  not  possible  to  explain  the 
radar  observations  in  the  direction  perpendicular  to  the  drift  motion  by  using 
the  above-mentioned  instability  mechanisms.  In  addition,  in  the  analysis  of 
those  instabilities,  energy  equations  are  not  incorporated  because  the  plasma 
temperature  perturbations  can  be  ignored  in  those  instability  processes. 

When  the  ohmic  dissipation  of  electrojet  currents  is  taken  into  account,  we 
find  that  a  new  instability  caused  by  the  thermal  effect  can  be  excited.  The 
physical  process  of  the  instability  is  through  the  modification  of  the  electron- 
neutral  collision  frequency  due  to  the  electron  temperature  perturbation  in  the 
electro  Jet.  The  collision  frequency  increases  with  the  temperature  perturbation. 
This,  in  turn,  further  increases  the  Joule  loss  os  the  electrojet  currents.  A  ther¬ 
mal  instability  is  thus  excited  through  such  a  positive  feedback  process.  A 
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dispersion  relation  of  the  instability  is  derived,  from  which  the  threshold  elec¬ 
tro  Jet  current  and  the  growth  rate  of  the  instability  are  determined.  It  is  shown 
that  the  analyzed  instability  mechanism  can  indeed  produce  those  as  observed 
relatively  large-scale  (greater  than  tens  of  meters)  E  region  irregular  structures 
which  are  polarized  perpendicular  to  the  drift  motion.  Details  of  this  work  are 
enclosed  in  the  Appendix. 

(D)  Overshoot  of  HF-Enhanced  Plasma  Lines  due  to  Resonance 
Broadening  Effects 

One  of  the  most  reproducible  phenomena  observed  in  the  Arecibo  ionos¬ 
pheric  heating  experiments  is  the  so-called  main  plasma  line  overshoot.  In  this 
work,  a  theoretical  model  is  developed  to  explain  the  observations.  This  model 
is  based  on  the  mode  suppression  process  introduced  by  the  resonance 
broadening  effect.  It  is  expected  that  the  presence  of  instabilities  gives  rise  to 
perturbations  on  the  phase  space  orbits  of  electrons.  The  result  of  the 
incoherent  scattering  of  electron  orbits  by  the  total  excited  waves  leads  to  elec¬ 
tron  diffusion  in  the  velocity  space  along  the  magnetic  field  together  with  the 
cross  field  diffusion  in  the  spatial  space.  These  diffusion  processes  thus  appear 
as  an  enhanced  viscosity  to  the  electron  motion  and  therefore  broaden  the  reso¬ 
nance  interaction  between  the  electrons  and  waves.  Consequently,  the  increase 
of  the  anomalous  damping  on  one  plasma  line  can  be  the  result  of  the  growth 
of  other  lines  in  the  same  region.  This  leads  to  the  mode-mode  competition 
and  hence,  the  spectral  lines  having  smaller  growth  rate  will  likely  be 
suppressed  by  the  presence  or  larger  growth  rate  lines.  Such  a  mode  competi¬ 
tion  process  leading  to  the  overshoot  of  HF-enhanced  plasma  lines  has  been 
described  in  terms  of  a  rate  equation,  which  determines  the  temporal  evolution 
of  the  parametrically  excited  plasma  lines.  The  result  of  the  modal  equation 
agrees  favorably  with  the  observations  in  the  Arecibo  heating  experiments. 
The  detail  of  this  work  is  also  enclosed  in  the  Appendix. 
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IV.  Appendix 

Reprints,  preprints  and  submitted  papers  attached  to  this  report  include: 

1.  Spectral  Characteristics  of  Hydromagnetic  Waves  in  the  Magnetosphere. 

2.  Thermal  Filamentation  Instability  Driven  by  the  Auroral  Electrojet 
Current. 

3.  Filamentation  Instability  of  Large  Amplitude  Alfven  Waves. 

4.  The  Stability  Analysis  of  a  Finte  Difference  Scheme  for  the  Coupled 
Hydromagnetic  Waves  Equations  in  the  Dipole  Model  of  the  Magneto¬ 
sphere. 

5.  A  New  Interpretation  of  Plasma-Line  Overshoot  Phenomena. 

6.  A  Theoretical  Model  of  Artificial  Spread  F  Echoes. 

7.  Simulanteous  Excitation  of  Large-Scale  Geomagnetic  Field  Fluctuations 
and  Plasma  Density  Irregularities  by  Powerful  Radio  Waves. 

8.  Resonant  Electron  Diffusion  as  a  Saturation  Process  of  the  Synchrotron 
Maser  Instability. 

9.  Parametric  Excitation  of  Whistler  Waves  by  HF  Heater. 

10.  On  the  Resonant  Ionospheric  Heating  at  the  Electron  Gyrofrequency 

11.  Enhanced  Ionospheric  Modifications  by  the  Combined  Operation  of  HF 
and  VLF  Heaters. 


J.  Plasma  Physics  (1987).  val  38.  part  2.  pp  235-243 
Printed  in  Great  Britain 


235 


Spectral  characteristics  of  hydromagnetic  waves 
in  the  magnetosphere 

By  S.  P.  KUO 

Polytechnic  University,  Route  1(0,  Farmingdale,  NY  11735 

M.  C.  LEE 

Massachusetts  Institute  of  Technology,  Cambridge,  Massachusetts  02139 

AND  A.  WOLFE 

New  York  City  Technical  College  of  the  City  University  of  New  York, 

Brooklyn,  NY  11201 

(Received  3  February  1987) 

This  work  is  intended  to  explain  why  the  resonant  response  of  the 
magnetosphere  prefers  to  have  discrete  frequencies.  Using  a  cylindrical  model 
for  the  outer  magnetosphere  with  a  plasma  density  profile  proportional  to 
1/r3,  we  show  that  the  eigenequation  characterizing  the  eigenmodes  of  the 
hydromagnetic  waves  in  this  model  has  two  turning  points  along  the  radial  axis. 
The  locations  of  the  turning  points  depend  upon  the  values  of  the  eigenperiod 
and  the  associated  east-west  wavenumber  of  the  eigenmode.  The  energy 
spectrum  of  the  excited  cavity  modes  is  seen  to  have  sharp  peaks  at  discrete 
frequencies  when  the  surface  perturbations  have  a  uniform  spectrum  in 
the  frequency  range  of  interest.  We,  therefore,  have  also  shown  that  only  the 
discrete  set  of  the  magnetospheric  cavity  eigenmodes  can  efficiently  couple  the 
perturbations  excited  on  the  boundary  of  the  magnetosphere  to  the  field-line 
resonant  mode  excited  inside  the  inner  turning  point  of  the  cavity  eigenmode. 
The  most  likely  values  of  east-west  wavenumbers  and  wave  period  range  are 
determined. 


1.  Introduction 

During  the  past  two  decades,  there  has  been  a  surge  of  interest  in 
understanding  geomagnetic  micropulsations.  Vigorous  observational  work 
(Allan  &  Poulter  1984;  Walker  &.  Greenwald  1981)  and  theoretical  research 
(Radoski  1970,  1971,  1973,  1974;  Chen  &  Hasegawa  1974;  Southwood  & 
Hughes  1983;  Rostoker  1979;  Lanzcrotti  &  Southwood  1979)  have  been  carried 
out  on  hydromagnetic  waves  and  geomagnetic  pulsations.  It  is  generally 
believed  that  micropulsations  arc  caused  by  hydromagnetic  waves  which  can  be 
generated  near  the  boundaries  of  the  magnetosphere,  internally  or  externally 
depending  on  local  time  and  on  the  frequency  ranges  of  the  pulsations.  For 
example,  hydromagnetic  waves  of  poloidal  and  toroidal  modes  transmitted 
along  geomagnetic  field  lines  can  be  excited  during  geomagnetic  storms  on  the 
nightside  and  they  can  be  excited  by  the  Kclvin-Helmholtz  instability  during 
quiet  times  on  the  dayside  magnetopause.  These  waves  then  give  rise  to 
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signatures  of  geomagnetic  pulsations  which  can  also  be  used  as  indicators  of  the 
changes  in  the  magnetosphcric  configuration.  The  resonant  field  line  picture 
(Chen  &  Hasegawa  1674:  Southwood  1974)  explains  successfully  the  variation 
of  the  amplitude  of  the  geomagnetic  pulsation  with  latitude.  This  model  also 
explains  very  well  the  associated  latitudinal  variation  of  the  wave  polarization 
(Samson,  Jackobs  &  Rostoker  1971;  Lanzerotti  el  at.  1974a, 6).  However,  the 
above  enumerated  work  does  not  explain  why  the  resonant  response  of  the 
magnetosphere  is  often  observed  to  have  a  discrete  spectrum  of  frequencies 
(Kokubum  &  Nagata  1965;  Stuart,  Sherwood  &  Macintosh  1971 ;  Takahashi  & 
McPherron  1982).  Moreover,  the  solution  of  the  field-line  resonance  mode 
predicts  spatially  dependent  eigenfrequencies  and  thus  contains  mathematical 
inconsistencies  in  the  analysis  (Radoski  &  McClay  1967). 

We  show  in  this  paper  that  the  cavity  modes  of  the  magnetosphere 
characterized  with  a  discrete  spectrum  of  frequencies  can  play  a  dominant  role 
in  geomagnetic  ULF  pulsations.  The  main  objection  to  the  concept  of  cavity 
modes  in  the  past  arises  from  the  lack  of  experimental  evidence  of  correlation 
between  the  frequency  spectra  of  the  dayside  and  nightside  signals.  This  fact  is 
inconsistent  with  the  assumption  of  the  standing  wave  pattern  around  the 
equator  (i.e.  along  the  east-west  direction)  for  a  spherical  cavity.  However,  this 
problem  can  be  resolved  if  one  takes  into  account  the  geometry  of  open  field 
lines  in  the  nightside  region,  namely,  assuming  that  the  open  field-line  region  is 
a  perfect  absorber  of  the  surface  perturbations  that  are  excited  in  the  dayside 
region  of  the  magnetopause.  Therefore,  the  wave  functions  in  the  east-west 
direction  should  have  a  travelling  wave  form  instead.  But  it  poses  another 
problem :  if  the  source  of  the  cavity  modes  has  a  broad  and  continuous 
wavenumber  spectrum  in  the  east-west  direction,  each  eigenmode  then  has  a 
finite  band-width  and  becomes  degenerate.  The  overlapping  of  the  eigen- 
frequency  spectrum  due  to  the  continuous  degeneracy  of  the  eigenmodes  will 
eventually  smear  out  the  discrete  nature  of  the  individual  cavity  mode. 

In  order  to  resolve  this  problem,  we  propose  a  model  that  considers  the 
magnetospheric  cavity  to  be  a  resonant  filter.  Then,  only  discrete  parts  of  a 
continuous  spectrum  of  the  perturbations  that  are  excited  in  the  magnetopause 
can  couple  efficiently  to  the  field-line  resonance  modes  located  inside  the 
magnetosphere  (Hasegawa,  Tsui  &  Assis  1983).  In  the  formulation  of  the 
theory,  we  use  a  cylindrical  model  of  the  magnetosphere  (Radoski  1970,  1976) 
to  analyse  the  eigenequation  of  the  resonant  ULF  waves.  We  will  first  show 
that  the  linearized  MHD  equations  can  be  combined  into  an  eigenequation 
lorm,  namely,  a  single  second-order  differential  equation  with  non-uniform 
coefficients.  When  a  plasma  mass  density  profile  proportional  to  1/r3 
(Cummings,  O’Sullivan  &  Coleman  1969;  Park,  Carpenter  &  Wiggin  1978)  is 
used,  it  is  seen  that  this  equation  has,  in  general,  two  real  turning  points  along 
the  positive  axis.  Therefore,  the  continuity  conditions  of  the  eigenfunctions  at 
the  turning  points  will  restrict  the  eigenvalues  of  the  equation  to  be  discrete  for 
a  fixed  cast-west  wavenumber  and  a  fixed  azimuthal  mode  number.  This  is 
essentially  the  same  concept  that  eigenvalues  become  discrete  in  a  potential 
well  as  used  in  quantum  mechanics.  The  cigcnmodc  solutions  are  undamped 
functions  inside  the  two  turning  points  in  the  collisionless  case.  Additional 
constraints  associated  with  the  preferential  cavity  modes  have  been  taken  into 
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Figure  1.  Cylindrical  magnetospheric  model  used  to  study  the  eigenfrequencies  of  the 
resonant  ULF  waves.  B  =  B(r)<}t.  The  z  direction,  normal  to  the  plane  of  the  figure,  is 
analogous  to  longitude. 

account  based  on  the  following  facts  :  (i)  the  locations  of  the  turning  points  are 
a  function  of  both  the  eigenvalues  and  the  corresponding  east-west  wave- 
numbers,  and  (ii)  the  perturbations  excited  in  the  magnetopause  can  be 
efficiently  coupled  to  the  inner  magnetosphere  only  when  the  outer  turning  point 
is  located  near  the  outer  boundary.  We  shall  show  in  the  following  analysis  that 
the  most  likely  values  of  east-west  wavenumbers  and  wave  period  range  can  be 
determined  with  the  present  model. 

2.  Theory  and  analysis 

In  a  cold  plasma,  the  hydromagnetic  wave  equations  derived  from  the 
Maxwell’s  equations  and  the  fluid  equations  have  the  general  form  of 

^  B°  x  {V  x  [V  x  (T)  x  B°)]}  =  0,  (1) 

where  tj  is  the  wave  plasma  displacement  vector,  B0  is  the  background 
geomagnetic  field,  and  £0  is  the  mass  density  of  the  plasma. 

This  wave  equation  will  be  solved  for  a  cylindrical  model  of  the 
magnetosphere.  In  this  model,  the  earth  becomes  a  flat  plane  and  the  magnetic 
field  lines  are  semicircles  as  shown  in  figure  1.  In  terms  of  the  cylindrical  co¬ 
ordinates  r,  $4,  z,  the  geomagnetic  field  is  expressed  as  B0  =  ifrBJr,  where  , B,  is 
a  constant.  It  is  found  that  (1)  has  only  two  non  trivial  components  (2  and  r) 
which  correspond  to  the  eigencquations  of  the  toroidal  mode  and  the  poloidal 
mode  of  the  hydromagnetic  waves,  respectively.  These  two  equations  will  be 
shown  to  couple  to  each  other.  If  we  assume  that  the  transverse  wave  electric 
fields  vanish  at  the  reflecting  points  on  the  conducting  earth  at  <f>  =  0,  n  it  would 
be  reasonable  to  use  a  standing  wave  expression  sinrt$4  for  the  azimuthal 
dependence  of  the  perturbations,  where  the  integer  n  measures  the  numbers  of 
half-waves  along  a  field  line.  The  2  direction  corresponds  to  the  longitude.  We 
assume  that  the  surface  perturbations  excited  in  the  daytime  regions  of  the 
magnetopause  disappear  when  they  propagate  toward  the  nightsidc  regions 
because  of  being  guided  by  open  field  lines.  A  total  absorption  boundary 
condition  then  requires  solutions  along  the  longitude  to  have  the  travelling 
wave  form  asexp[»(I,2  — w()]  with  continuous  kt.  Thus  we  may  assume  that  the 
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wave  perturbations  take  the  form  T]r  =  [F(r)/.R0(r)]  siting  exp  [i{ktz  —  wt)\  and 
i]z  =  [H(r)/ Bv(r)]  sin rup  exp[t(fczz  —  wt)],  and  the  two  coupled  eigenmode  equa¬ 
tions  have  the  following  forms  (Radoski  1970,  1976): 


< oY  n 2  Id  d 
vAJ  r2  +  r  dr  r  dr 


F(f)  =  —  ikz  —  -y-  rH (r), 
r  ar 


(2) 


W(r)  = 


-t*2-F(r), 


(3) 


where  =  R0/(47r£0)i  is  the  Alfven  speed,  and  H(r)  and  F(r)  are  functions 
characterizing  the  radial  dependence  of  the  amplitudes  of  the  torodiai  and 
poloidal  modes,  respectively. 

From  (2)  and  (3),  an  eigenequation  for  the  shear  magnetic  perturbation 
G(r)  =  ilczrH(r)  +  r(d/dr)F(r)  is  derived  to  be 


d2 

dr2 


G— 


^  In  (rfc2)  ^-G  +  (A*-Jfc*)0  ==0, 


(4) 


where  the  notation  A;2  =  ((i)/vA )2  —  n2/r2  is  used. 

Letting  (7(r)  =  (rA:2)1 <?,(»•),  equation  (4)  then  reduces  to  a  wave  equation  with 
a  spatially  dependent  wave  vector  given  by 

£iGl(r)+f<r)Gl(r)=  0,  (5) 


where  /(r)  =  A-2  —  A:2  —  J  a2  +  ^  ^  a 


and 


a  =  —  In  (rA-2). 
dr 


We  now  model  the  outer  magnetosphere  with  a  plasma  density  profile 
£o  —  h/**  f°r  ^  r  ^  UR,,  where  is  a  constant  and  Re  is  the  earthradius. 
In  the  present  paper,  the  sharp  plasma  density  gradient  located  at  the 
plasmapause  (r  %  2  — 3R,)  is  not  taken  into  account.  Thus 


and 


/(»•>  = 


2n*  +  sr  t  _  n2  1 
4(n2  — sr)2  ‘  r2  +  4r2’ 


(6) 


where  s  =  47t£,  w2//i2  is  a  parameter  proportional  to  the  squared  eigen- 
frequency  w. 

As  shown  in  figure  2,  /(r)  has  two  zeros  located  at  two  points  r,  and  rf  on 
the  positive  real  axis.  It  implies  that  (5)  has  two  real  positive  turning  points 
r,  and  r,.  At  the  turning  points,  the  wavenumber  /(r)i  becomes  zero,  thus  each 
of  them  defines  a  boundary  separating  the  oscillatory  and  non-osciliatorv  parts 
of  the  solution  of  the  wave  equation.  In  other  words,  the  wavefunction  begins 
to  change  its  characteristics  at  the  turning  point.  However,  this  change  must  bo 
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Figure  2.  Spectral  variations  of  the  square  of  the  wavenumber /(r)  from  r  =  2  Re  to 
r  —  HR'  where  n  =  1,  lct  =  0  32 /R,  and  e  =  1-12//?,  are  used. 

compromised  by  the  continuity  conditions  of  the  wavefunction  and  its  first 
derivative  at  each  of  the  turning  points  so  that  the  second  derivative  of  the 
wavefunction  can  remain  finite  as  imposed  by  the  considered  wave  equation. 
The  position  of  the  turning  points  r,  and  r2  depends  on  the  values  of  the 
parameters  s,  n  and  kt.  Furthermore  f(r)  is  negative  in  the  region  outside  the 
turning  points.  We  find  that  (5)  has  solutions  for  G(r)  which  are  well  behaved, 
namely,  they  decay  exponentially  away  from  the  turning  points  in  that  region. 
Thus  spatially  localized  modes  exist  provided  an  eigenvalue  s  (i.e.  cu)  can  be 
found  such  that  the  eigenvalue  equation  for  s  (or  w),  defined  by  the  Bohr- 
Sommerfeld  quantization  condition  (Landan  &  Lifshitz  1965) 

(m  +  |)n  «  £  fhr)dr,  (7) 

can  be  satisfied,  where  m  =  0, 1,2,...  stands  for  the  number  of  nodes  of  the 
eigenfunction  G,(r)  within  the  two  turning  points  r,  and  rt. 

Equation  (7)  is  solved  numerically  for  s.  Since  r,  and  rt  are  also  functions  of  s, 
equation  (7)  has  to  be  solved  self-consistently.  Presented  in  figure  3(a)  and  (6) 
are  examples  of  spatially  localized  eigenfunctions  determined  by  (5)  where  (7) 
has  been  used  to  determine  first  the  approximate  eigenvalues.  For  the  given 
values  of  n  and  kt,  we  find  that  (7)  only  has  discrete  solutions.  Moreover,  the  two 
turning  points  have  to  be  located  inside  the  two  boundaries  taken  at  2 /?,  and 
11/?,  of  the  modelled  magnetosphere  in  order  to  sustain  spatially  localized 
modes.  This,  in  turn,  determines  the  lower  bound  ktl  and  the  upper  bound  klt 
of  the  wavenumber  kt.  In  other  words,  for  a  given  integral  value  of  n,kt  must 
be  within  the  interval  (ktl,kti)  so  that  the  turning  points  of  (5)  determined  by 
each  of  the  solutions  of  (7)  satisfy  the  conditions  r,  >  2/?,  and  r,  ^  11  Re. 
Figure  4  displays  the  plots  of  /(r)  for  the  tw  o  extreme  cases  kt  =  ktl  and 
kt  =  kt j,  respectively.  We  observe  that  r,  >  2 Re  and  r,  =  11/?,  in  the  first  case 
of  Jfc,  =  ktl  and  s  —  s,,  and  r,  =  2/?,  and  r,  <  1 1/?,  in  the  second  case  of  kt  =  kn 
and  s  =  Sj.  If  the  resonant  ULF  waves  in  the  magnetosphere  are  assumed  to 
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n 

m 

7,  (sec) 

7 j  (sec 

1 

0 

217 

130  3 

18  2 

119  5 

2 

0 

15  1 

872 

90 

67  1 

3 

0 

119 

662 

6-0 

470 

4 

0 

97 

53  1 

4  4 

35-8 

4 

1 

71 

40-9 

60 

37-7 

5 

0 

8-3 

44  7 

35 

29  1 

5 

1 

6-2 

354 

4-6 

306 

6 

0 

7-3 

38-8 

2-8 

24  1 

6 

1 

56 

314 

36 

253 

Table  1.  Upper  and  lower  bounds  of  the  wave  periods  and  east-west  wavelengths  of  the 
eigenmodes  for  azimuthal  mode  number  n  from  1  to  6. 


rad  s"1 

Fig L re  5.  Energy  spectrum  of  the  magnetospheric  cavity  modes  excited  by  a  source  of 
perturbation  with  a  uniform  spectrum  in  the  magnetopause. 


arise  from  the  perturbations  in  the  magnetopause,  i.e.  excited  by  the  Kelvin- 
Helmholtz  instability  (Chen  &  Hasegawa  1974 ;  Southwood  1974),  it  is  quite 
obvious  that  the  perturbations  with  parameters  corresponding  to  case  one  can 
couple  into  the  magnetosphere  more  efficiently  than  that  corresponding  to  case 
two.  This  is  because  in  case  one  the  outer  turning  point  is  right  in  the 
magnetopause,  hence  the  perturbation  can  be  directly  propagated  into  the 
magnetosphere  without  tunnelling  through  a  barrier  as  in  case  two.  Table  1 
lists  the  computed  eigenperiod  bands  for  the  azimuthal  mode  number  n 
from  1  to  6.  The  following  L  =  4  reference  parameters:  B0  =  4  86  x  10-’  G,  £0  = 
6  68  x  10~”  g  cm-*,  Rt  =  6  4  x  10*  cm  are  used  in  the  computation  (Helliwell  & 
Inan  1982).  It  shows  that  the  eigenperiod  bands  for  n  =  3  to  6  overlap  with  one 
another.  This  is  because  each  (n,  m)  mode  has  a  finite  fc,  band.  Each  eigenmode 
then  has  a  finite  band-width  and  becomes  degenerate  in  the  overlapped  range 
of  the  wave  spectrum.  However,  these  modes  do  not  have  the  same  intensities 
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because  their  turning  points  are  different.  Those  having  turning  points  away 
from  the  magnetopause,  that  is  the  source  location  of  perturbations,  have 
weaker  initial  intensities  weighted  by  an  exponentially  decaying  factor, 

expj-j" 

In  this  factor,  /  is  defined  by  (6),  11  Re  is  the  radial  distance  of  the 
magnetopause  from  the  centre  of  the  earth  and  r,  represents  the  radial  distance 
of  the  turning  points  concerned.  As  illustrated  in  figure  5  following  table  1 ,  even 
when  the  source  perturbations  in  the  magnetopause  have  a  uniform  spectrum 
(expressed  in  an  arbitrary  unit)  in  the  frequency  range  of  interest,  the  energy 
spectrum  of  the  excited  cavity  modes  can  have  sharp  peaks  at  discrete 
frequencies  in  agreement  with  observations.  The  discrete  nature  of  the  cavity 
mode  spectrum  can  be,  therefore,  viewed  as  the  consequence  of  the  following 
fact.  A  potential  barrier  is  imposed  on  the  coupling  of  magnetospheric 
perturbations  (at  the  magnetopause)  into  the  cavity  modes  as  the  outer  turning 
points  of  the  cavity  modes  do  not  match  the  location  of  the  magnetopause. 
Therefore,  in  table  1,  the  values  ofrt  and  An  are  considered  to  be  the  most  likely 
periods  and  east-west  wavelengths  of  the  field-line  resonances  excited  in  the 
magnetosphere. 


3.  Summary 

We  have  theoretically  explained  the  discrete  nature  of  ULF  pulsations 
observed  in  the  magnetosphere  by  using  a  cylindrical  model  (Radoski  1970, 
1976).  A  second-order  wave  equation  is  derived  to  characterize  the  coupling 
between  the  poloidal  and  toroidal  modes  of  the  hydromagnetic  waves.  When 
the  experimentally  determined  plasma  density  radial  profile  of  the  ‘  1/r3  ’  form 
(Cummings  el  al.  1969;  Park  el  al.  1978)  is  used,  it  is  found  that  the  wave 
equation  has  two  turning  points  along  the  radial  axis.  Undamped  eigenfunctions 
with  discrete  eigenvalues  are  found  in  our  collisionless  case.  They  are 
conceptually  analogous  to  the  case  that  has  discrete  eigenstates  when  a  particle 
is  trapped  by  a  potential  well  as  discussed  in  quantum  mechanics.  However, 
there  is  a  difference  between  the  two  cases.  In  the  present  case,  f(r)  cannot 
simply  be  expressed  by  the  difference  of  a  constant  eigenenergy  and  a  fixed 
potential  function  as  in  quantum  mechanism.  This  is  because  the  eigenvalues 
(or  <u)  are  multiplied  by  a  spatially  dependent  function.  Therefore,  the  potential 
well  determined  by /(r)  varies  with  each  discrete  eigenstate.  The  eigenfunctions 
as  shown  in  figure  3(a)  and  (6)  decay  rapidly  toward  the  ‘equator’.  However, 
the  eigenfunctions  still  remain  finite  at  r  =  1-6 R,  in  the  case  demonstrated  in 
figure  3(6).  This  means  that  the  pulsations  can  still  be  observed  at  the  surface 
of  the  earth  in  the  polar  regions  where  they  arc  a  distance  {wfl,  away  from  the 
equator.  We  have  calculated  the  eigenperiods  r,  of  the  most  readily  excited 
modes  and  the  most  likely  values  of  the  east-west  wavenumbers  Atl  as 
tabulated  in  table  1  for  illustrative  purposes.  The  wave  eigenperiods  are 
calculated  to  lie  in  the  Pc  3-4  range. 

Further  work  should  include  the  sharp  plasma  mass  density  gradient 
occurring  in  the  inner  magnetosphere  (plasmapause)  together  with  a  realistic 
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plasma  variation  for  the  outer  magnetosphere.  These  results  using  a  cylindrical 
model  should  then  be  modified  with  a  more  realistic  dipole  field  model  of  the 
earth.  Transient  time  damping  of  the  cavity  modes  by  hot  plasmas  due  to 
kinetic  effects  should  also  be  included  phenomenologically  in  the  wave  equation. 
In  this  case,  the  cavity  modes  having  complex  frequencies  become  damped 
quasi-modes. 
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A  thermal  instability  leading  to  the  filamentation  of  auroral  electro-jet  currents  and  giving  nse  to 
purely  growing  magnetic  Held-aligned  density  irregularities  in  the  £  region  of  the  high-latitude  iono¬ 
sphere  is  investigated.  The  physical  process  of  the  instability  is  through  the  modification  of  the  electron- 
neutral  collision  frequency  due  to  the  electron  temperature  perturbation  in  the  electrojet  A  dispersion 
relation  of  the  instability  is  derived,  from  which  the  threshold  electrojet  current  and  the  growth  rate  of 
the  instability  are  determined.  It  is  found  that  they  become  independent  of  the  scale  sizes  of  the  irregu¬ 
larities  for  scale  sizes  larger  than  about  13  m.  The  proposed  instability  can,  thus,  be  considered  to  be  one 
of  the  mechanisms  responsible  for  those  observed  relatively  large-scale  £  region  irregularities  (Pfaff  et  al., 
1984). 


1.  Introduction 

The  high-latitude  ionosphere  is  constantly  structured  in  the 
ionospheric  irregularities,  (see,  for  example,  Kelley  and  Mozer 
[1973],  Phelps  and  Sagalyn  [1976],  Kino  et  at.  [1978],  Vickrey 
et  al.  [1980],  Weber  et  al.  [1984],  and  Weber  and  Buchan 
[1985].  The  appearance  of  these  density  irregularities  has  gen¬ 
erally  been  attributed  to  the  excitation  of  plasma  instabilities. 
These  include  the  Farley-Buneman  two  stream  instability 
[Farley,  1963],  and  the  E  x  B  and  current-convective  instabil¬ 
ities  [Ossakow  and  Chaturvedi,  1979;  Keskinen  and  Ossakow, 
1983;  Mitchell  et  al.,  1985].  Recently,  Weber  et  al.  [1984]  and 
Cerisier  et  al.  [1985]  have  invoked  the  E  x  B  instability,  and 
Das  and  Das  [1983]  and  Lee  [1984]  have  proposed  thermal 
instability  to  explain  large  scale  density  fluctuations  in  the  F 
layer  of  the  high  latitude  ionosphere,  and  Chaturvedi  et  al. 
[1987]  have  considered  the  parallel  current  effects  on  the 
onset  of  the  Farley-Buneman  instability  in  the  auroral  E 
region.  A  comprehensive  review  of  ionospheric  irregularities 
excited  by  some  instabilities  has  been  published  in  Review  of 
Geophysics  by  Fejer  and  Kelley  [1980],  and  in  Radio  Science 
by  Keskinen  and  Ossakow  [1983].  Both  large  and  small-scale 
irregularities  have  been  observed  from  the  £  and  F  layers  of 
the  high-latitude  ionosphere.  The  Farley-Buneman  instability 
driven  by  electrojet  currents  is  known  to  be  effective  in  gener¬ 
ating  meter  and  shorter-scale  irregularities  of  type  /  in  the  £ 
region.  This  instability  has  been  extensively  investigated  by 
radar  techniques.  The  irregularities  of  type  II  are  of  slightly 
larger  scale  sizes  and  are  generated  through  the  E  x  B  gradi¬ 
ent  drift  instability.  However,  the  growth  rate  of  the  two 
stream  instability  as  well  as  of  the  gradient  drift  instability  is  a 
function  of  the  angle  between  the  wavevector  and  the  direc¬ 
tion  of  the  electron  drift.  The  growth  rate  becomes  zero  when 
the  wavevector  is  perpendicular  to  the  drift  direction.  It  is 
therefore  not  able  to  explain  the  radar  observations  in  the 
direction  perpendicular  to  the  drift  motion  by  using  the  linear 
two  stream  instability  mechanism.  A  two-dimensional  nonlin- 
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ear  theory  [Sudan  et  al.,  1973]  has,  thus,  been  proposed  to 
explain  the  presence  of  the  type  II  gradient  drift  irregularities 
with  wavelengths  and  propagation  directions  that  are  predic¬ 
ted  to  be  linearly  damped.  Furthermore,  the  temporal  and 
spatial  power  spectrum  of  the  irregularities  in  the  equatorial 
electrojet  have  also  been  computed  [Keskinen  et  al.,  1979] 
from  simulations  of  the  two-dimensional  model.  The  results 
are  found  to  compare  favorably  with  several  9-m  type  II  radar 
spectra  taken  in  the  daytime  equatorial  electrojet  [Balsley  and 
Farley,  1971],  However,  in  the  analyses  of  these  instabilities, 
energy  equations  are  not  incorporated  because  plasma  tem¬ 
perature  perturbations  can  be  ignored  in  these  instability  pro¬ 
cesses. 

If  the  ohmic  dissipation  of  electrojet  currents  is  taken  into 
account,  we  expect  that  a  new  instability  caused  by  the  ther¬ 
mal  effect  may  be  excited.  In  this  paper,  we  study  this  new 
mechanism  for  the  generation  of  irregularities  in  the  auroral 
electrojet  under  the  condition  that  the  electron  drift  speed 
exceeds  3  times  the  ion  acoustic  speed.  This  thermal  insta¬ 
bility,  will  be  shown  also  to  produce  relatively  large-scale 
(greater  than  tens  of  meters)  £  region  irregular  structures 
which  are  polarized  perpendicular  to  the  drift  motion.  It  is 
then  contributed  additively  to  the  nonlinear  gradient  drift  in¬ 
stability.  This  instability  is  first  analyzed  analytically  and  de¬ 
scribed  in  section  2,  in  which  the  dispersion  relation  of  the 
instability  is  presented.  This  dispersion  relation  is  analyzed  to 
determine  the  threshold  and  growth  rate  of  the  instability.  The 
results  are  presented  in  section  3.  Conclusions  and  general 
remarks  are  given  in  section  4. 

2.  Dispersion  Relation 

We  study  the  thermal  instability  excited  by  the  background 
efectrojet  current.  The  vector  orientation  of  Figure  1  is  con¬ 
sidered.  where  E0  =  v£0  is  the  dc  driving  field  of  the  electrojet 
current  and  B0  =  —:Bn  is  the  geomagnetic  field.  Hence  the 
background  Pederson  <x  component)  and  Hall  (y  component) 
current  densities  are  given  respectively  by 

fro  =-  -  '\o«>  “  «") 

[v,0,2  +  n,2)  +  (m  +  O,2)]£0  (I) 
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F  ig.  I.  Kilamentation  of  electrojet  current  by  thermal  instability. 


Jy0  =  -‘'«o(<',oy  -  ‘ioy>  =  -(nae2im) 

[0,/(v/  +  a/)  -  ( m/MfQ'Av (2  +  n,2)]£0  (2) 

where  m( XI v^v,)  and  Q,(fl,)  are  the  electron  (ion)  mass,  un¬ 
perturbed  collision  frequency  and  cyclotron  frequency  respec¬ 
tively;  n„  is  the  background  plasma  density;  v,0  = 
-  [(c£0/m)/(v,2  +  n,2)](xv,  -  yft,)  and  vl0  ==  [(e£0/M)/(v,2 
+  fti2)]|xvi  +  >«,). 

In  the  presence  of  thermal  fluctuations,  the  plasma  density 
and  the  electron  temperature  of  the  electrojet  are  perturbed. 
In  terms  of  the  unperturbed  quantities  n0  and  Te0 ,  the  total 
density  n  and  the  electron  temperature  Tr  are  expressed  as 
n  =  na  +  tin  and  Te  =  Tr„  +  <57,,  where  the  quasineutrality  is 
assumed  so  that  <5 n,  -  Sn  =  Sn ,  and  Sn  and  <57,  are  the  density 
and  temperature  pertubations,  respectively.  These  pertur¬ 
bations  are,  in  turn,  to  introduce  modifications  on  electron 
collision  frequency  and  the  plasma  conductivity.  Since  the 
electron  collisions  are  dominated  by  the  electron  neutral  col¬ 
lision  in  the  £  region,  the  collision  frequency  v,  is  proportional 
to  the  square  root  of  the  electron  temperature  and  can  be 
expanded  up  to  the  lowest  order  temperature  perturbation  as 
v,  ^  v,(l  +  $STr/Te0).  The  conductivity  will  also  be  expanded 
up  to  the  lowest  order  density  and  temperature  perturbations 
in  a  similar  way.  With  the  aid  of  these  expansions  and  as¬ 
suming  that  the  linearized  perturbations  bear  the  general  form 
SP  =  SPe* Vyl  having  real  positive  y  as  the  growth  rate  of  the 
field-aligned  purely  growing  mode  and  k  =  xkx  +  yky  perpen¬ 
dicular  to  the  magnetic  field,  the  transport  equations  given  by 
Braginskii  [1965]  are  linearized  to  be 

k  ■  (<5v,  -  <5v,)  =  —  (<5n/n0)k  u0  (3) 

k  •  <5v,  =  (/>■  -  k  •  vj0X<5n/n0)  (4) 

obtained  from  the  electron  and  ion  continuity  equations, 
where  u0  =  v,0  -  vl0, 

v,„k  ■  <5v,  -  Q,k  x  <5v,  •  r  =  -  i(k2/mK<5T, 

+  Tef)Sn/n0)  -  jv,„k-  y„{STjTe0) 
n,k  •  <5v,  +  v,„k  x  <5v,  •  f  =  -  K„(k  *  V  ’  Te0)  (5) 


v,„k  •  (<5v,  -  <5v,j  +  (M/m)v,„k  ■6\l  -  0,k  x  (<5v,  -  <5v,)  •  £ 

=  -Hk2/m)[5T'+{Tt0  +  7]X<5n/n0)]-*v,.(k  ■  y^STJT^)  (6) 

ft,k  •  (<5v,  -  <5v f)  +  v,„k  x  (<5v,  -  ^Vj)  •  £  +  ( M/m)v ^k  x  <5 v  * 

=  ~  2  VJ£  •  k  x  (2) 

obtained  from  the  electron  and  ion  momentum  equations, 
where  =  v„  +  (m/M)v„  +  y  +  ik  •  vi0,  and 

(V  +  «k  •  T'0 k  • 

=  -v„(*  V/n.W,  -  (m/M.)v„(3  -  TJTr0)STc 

+  (2/3no0(O/,)  -  2 (m/MJyJT^  -  TJ6n/n0)  (8) 

derived  by  the  electron  thermal  diffusion  equation,  where 
vu  =  (Tr0/m)'12  and  the  differential  ohmic  heating  source 
S(Qfy)  =  S{J1/nal*\  with  <r1*  =  (1/<T  +  Mtl(2/ne2vJ' 1  and 
a  =  ne2/mvn  [ Braginskii ,  1965],  is  given  by 

—  ^Q/,)  =  («o/l’<.)2[2  v'.6T.- T'O(v„  +  {lfll/vlm)(6n/n0)'} 
no 

+  (m/n02e2Xv„  +  ClJClJvJSJ2  (9) 

We  can  further  express 

SJ2  =  2J02(6n/n0)  -  2en0J0  ■  (<5v,  -  <5v,) 

=  2J02(Sn/n0)  +  2e2n02{(k  •  u0)[k  •  (<5v,  -  <5u,)] 

+  (z  -k  x  u0)[z'  •  k  x  (<5v,  -  <5v,)]}/k2  (10) 

Solving  (5H7)  yields 

k  •  (<5v,  -  <5v()  =  {[(mv„  -  Mvjn,  +  mfi,c,](<5r^T,0) 

+  -  JWvJb,  -  mfl ,d,] 

•  (<5n/n0)}/[mv„  -  M v  J2  +  m2ft, 2]  (11) 

and 

k  x  (<5v,  -  <5vf)  •£  =  {[(mv,„  -  Mvjc,  -  m£ltaJ(STJTr0) 

-  [<«v«,  -  Mvill)di  +  mfi,b,] 
(«5n/n0)}/[(mv„  -  Mu,.)2  +  m2ft,2]  (12) 

where 

a,  =  -tk2T^-i»nv„k  •  v,0  +  [(Mv,n/m(v„2 +  fl,2)](av„ +  cQ,) 
h,  =  —ik2(Te0  +  7;)  +  [Mv(ivjm(v„2  +  n,2)]/> 

(13) 

<  i  =  -  •  k  *  >,o  +  tWvln/ni(v„2  +  Q,2)](cv„-an,) 

d,  =  [Mv,A/m(v„2  +  ft,2)]h 

and 

a  =  ik2Tr0  +  imv,nk  •  v,0 
b  =  ik2Tr0 
<'  =  i'f'V.k  x  o  •  ; 

From  (3)  and  (1 1),  we  have 

|k  •  u0  +  [(mi„  -  Mijh,  -  mQrd , ]/[(mv„  —  Mi’,,)2  +  m2I2, 2] } 
(<5n  n0)  =  -  ![mv„  -  Mv„)a,  +  mn,c,]/[(mv„  -  MvJ2 
+  m2n,2]H<5T,  ,Tr0)  (15) 

substituting  (9)  into  (8),  and  with  the  aid  of  (4),  (5).  (10),  (11), 
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(12),  and  (IS),  the  dispersion  relation  for  the  filamentation 
instability  of  the  electrojet  is  derived  to  be 

y  +  ik  •  v,o  -  i v„[(u0/t„)2  -  5 k2v,t2/nt2  -  3(m/M,M3  -  TJTt0)] 

-i  ■  f,o  +  k  *  *,o  ■’! 

=  f  v„{[«c2p„2/2t2,2)  +  (3m/2M„)(l  -  TJT^  +  d 
•  [(It  •  Uo /kvj2  -  i(u0/vuf]}l(mV'n  -  Mv„)ii,  + 

+  (1  +  OA/v„vJ[(z  •  k  x  u0)/k2vu2] 

{(It  •  u0)[(mv„  -  M vlm)c ,  - 

+  (b,c,  +  <J,d,)]}/{(k  •  u0)[(mv„  -  MvJ2  +  m2CJ,2] 

+  [('”*,„  ~  Mv>,  -  (16) 

This  equation  will  be  solved  in  the  following  to  determine 
the  threshold  electrojet  current  and  the  growth  rate  of  the 
instability,  and  the  filamentation  direction  of  the  electrojet. 
The  thermal  source  of  the  instability  is  given  by  the  ohmic  loss 
of  the  electrojet  current  manifested  by  terms  proportional  to 
vjt',0  —  t’,0r.  Collision  and  diffusion  losses  impose  a  thresh¬ 
old  condition  of  the  instability. 

3.  Results 

Equation  (16)  is  complex  which  can  then  be  broken  up  into 
two  real  equations.  While  the  real  part  of  (16)  determines  the 
growth  rate,  the  imaginary  part  of  (16)  defines  the  direction  of 
filamentation  instability.  These  two  equations  are  given,  re¬ 
spectively,  to  be 

y  =  K„[(u0A’J2  -  S*VA*  -  WMJI  -  TjiTrS 

+  3v»2vjll2/Q,){  -(*2t>„2/nr2) 

-  (3m/M„Xl  -  TJT'q)  +  (u0/OJ 
+  2(fl,ni2/vMvi„2)[(Ol/vl„Xk  •  n0)2 

—  (k  •  u0Xz  •  k  x  u0)]/k2t’„2} 

■  (f  •  k  x  v,0>/C*'m 2  +  2*k  ‘  uo>]  (17) 

and 

<k  •  *,<>)[ 1  “  3<vw/n,)2]  -  3<'V'n«J<z  •  k  x  v,0) 

s  i[v„2v(.2/nA<v,.2  +  a,2)] 

(1  +  OA/v„vJ(z'  •  k  x  v<-oKz  •  k  x  u0)/(k  •  u0) 

+  }(v JOX 1  +  *Wv„vJ(z  •  k  x  u0)  (18) 

where  we  have  approximated  v,„  by  vto,  i.e.,  |y  +  (m/M)vfll  +  ik 
•  \i0 1  «  v,„  is  assumed,  and  neglected  the  ion  thermal  effect, 
simplify  the  resultant  equations. 

In  the  quantitative  analysis,  we  generally  assume  uneven 
electron  temperature  (T,0)  and  neutral  temperature  (T„  ~- 
7j0  -  300  K).  and  the  following  E  region  parameters:  v,„  =  5 
x  I04s  v,.  =  3  x  102  s  ',  =  159  s  ’,  Qe  =  8.8  x  10‘ 

s’1.  M,(NO*)/m  =  5.52  x  10*,  and  r„  =  (Tf0/TJll26J4  x 
10*  m/s  will  be  used. 

We  first  analyze  (18)  to  determine  the  filamentation  direc¬ 
tion.  i.e.,  the  direction  of  k.  With  the  aid  of  (10)  and  (2),  and 
using  the  £  region  parameters  for  simplifying  the  result,  we 
obtained 

0  =  tan  Hk„kx)  *  tan-  ,[4v1>Qi/(3vi(IJ  -  fl,2)]  (19) 

From  (1)  and  (2).  the  direction  of  the  electrojet  current  is 


also  obtained  to  be 

0o  =  tan -  ‘(7vo/7,ot  -  -  tan '  '(vJO,)  (20) 

consequently,  it  is  found 

0  -  0O  £  -tan  '(3v/fl,)  -n/2  (21) 

We  have  thus  shown  that  the  result  of  thermal  instability  is 
developing  the  current  filaments  along  the  electrojet  current. 
The  excited  instability  is  both  magnetic  field  aligned  and  elec¬ 
trojet  current  aligned.  This  feature  is  the  characteristic  nature 
of  the  filamentation  instability  that  is  associated  with  purely 
growing  (i.e.,  zero  real  frequency)  mode.  If  the  excited  mode  is 
not  electrojet  current  aligned,  this  mode  ihen  has  finite  fre¬ 
quency  determined  by  k  •  u0.  It  is  found  that  higher  threshold 
is  required  for  excited  nonpurely  growing  mode.  Therefore,  the 
physical  reason  for  proposed  instability  to  be  electrojet 
current-aligned  in  addition  to  magnetic  field-aligned  is  that  it 
meets  the  minimim  threshold  for  its  excitation. 

We  next  use  (19)  for  the  direction  of  k,  (17)  then  becomes 

V  =  iv,„{(uo/eIf)2r(40O*  +  7Q(2vi.2  +  9v1/)/(flj2  +  v(„2) 
(16Q,2  +  90  +  IwS/Clfl/i  fl,2  +  v(„2)] 

-  5(ke,yn,)2  -  9(m/MJl  -  TJ3TJ 
[1  +  v„vi,3/n,n((vj,2  +  oi2)]}  (22) 

Using  the  E  region  parameters,  we  can  express  (22)  as 
V-  I.19[(u0/C,)2  — (Tf0/r.X162/A2)  — 2.52(1  -TJIT^  (23) 

where  c,  =  ~  500  (Tr0/T,)'<2  m/s. 

The  threshold  condition  obtained  from  (23)  by  letting  y  =  0 
is  found  to  be 

'  [(T,0/T„X  162/a2)  +  2.52(1  -  TJ3 T^v,  (24) 

It  is  a  functional  dependence  of  T^/T,.  In  observations, 
T»o/T„  is,  usually,  less  than  5  (R.  Pfaff,  private  communication, 
1987).  We  therefore  limit  our  analyses  in  the  following  for 
T'o/T.  <  5. 

Figure  2  presents  the  relation  between  the  threshold  drift 
speed  u0lh  and  the  scale  length  A  of  the  thermal  instability.  It 
shows  that  u0tJv,  quickly  approaches  to  a  minimum  about 
1.55  for  A  >  30  m  in  both  T^/T,  =  2  and  5  cases. 

The  dependences  of  the  growth  rate  y  on  the  scale  length 
are  shown  in  Figure  3  for  two  cases  in  corresponding  with 
Figure  2.  We  first  consider  the  case  T^/T,  —  2  and  assume 


Fig.  2.  Threshold  drift  speed  versus  the  scale  length  of  the  thermal 
instability,  where  Tr„,Tx  —  2  and  r,  -  707  m  s  (solid  line)  and 
Tr0  T„  =  5  and  r,  —  11 20  m  s  (dashed  line)  are  assumed. 
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Fig.  3.  Growth  rate  versus  the  scale  length  of  the  irregularity, 
where  the  solid  line  represents  the  case  £0  =  75  mV/m,  u0  =  1.5  km/s, 
and  T^'T  -  2,  and  the  dashed  line  stands  for  the  case  £0  =  150 
mV/m,  u0  -  3  km/s,  and  T^JT,  -  5. 


that  u0  =  1.5  km/s  corresponding  to  a  dc  driving  field  £0  =  75 
mV/m.  In  this  case,  only  perturbations  with  scale  lengths 
larger  than  11.5  m  can  be  excited.  The  growth  rate  quickly 
reaches  a  maximum  2.86  s”  1  for  2  >  30  m.  We  next  consider  a 
strong  electrojet  current  case  with  £0  =  150  mV/m  and  u0  =  3 
km/s.  strong  ohmic  heating  leads  to  large  T„IT,  and  a  ratio  5 
is  thus  assumed.  The  minimum  scale  length  of  the  instability 
increases  only  slightly  to  13  m.  However,  the  maximum 
growth  rate  increases  to  about  twice  the  previous  case. 

4.  Discussion  and  Conclusion 

We  have  shown  that  a  thermal  instability  can  be  excited 
when  the  electrojet  currents  are  intense  enough.  This  insta¬ 
bility  can  cause  the  filamentation  of  the  electrojet  currents  and 
gives  rise  to  tens  to  hundreds  of  meters  £  region  irregularities. 
The  physical  process  of  the  instability  is  through  the  modifi¬ 
cation  of  the  electron-neutral  collision  frequency  due  to  the 
electron  temperature  perturbation  in  the  electrojet.  The  col¬ 
lision  frequency  increases  with  the  temperature  perturbation. 
This,  in  turn,  further  increases  the  ohmic  loss  of  the  electrojet 
currents.  A  thermal  instability  is  thus  excited  through  such  a 
positive  feedback  process.  Here  we  have  considered  a  mag¬ 
netoplasma,  and  hence  the  heat  conduction  in  the  direction 
perpendicular  to  the  magnetic  field  is  greatly  reduced.  In  this 
case  instability  becomes  favorable  to  be  excited  in  magnetic 
field-align  structure  and  the  heat  transport  processes  are  es¬ 
sentially  determined  by  the  collisional  coupling  between  elec¬ 
trons  and  neutrals.  Electron  gas  absorbs  the  energy  from  the 
heat  source  and  dumps  a  fraction  of  it  to  the  neutrals.  For 
scale  sizes  larger  than  about  13  m,  the  diffusion  loss  can  be 
neglected  and  thus  the  threshold  and  the  growth  rate  of  the 
instability  become  independent  of  the  scale  sizes  of  the  irregu¬ 
larities. 

We  also  note  that  unequal  electron  and  neutral  temperature 
in  the  high  latitude  ionosphere  were  observed  [ Schlege I  and 
St. -Maurice.  1981;  Wickwar  et  at.  1981;  Stauning.  1984], 
Ohmic  dissipation  of  the  electrojet  current  may  be  partly  re¬ 
sponsible  for  the  elevated  electron  temperature.  High  thresh¬ 
olds  are  required  for  the  instability  because  of  the  additional 
damping  imposed  by  Te0  >  T„  condition  as  shown  in  Figure  2. 
The  threshold  can  be  exceeded  at  high  latitude  when  intense 
electrojet  currents  are  present.  Hundred-meter  scale  irregu¬ 
larities  can  be  excited  by  the  thermal  instability  in  less  than  a 
second.  Once  the  threshold  conditions  are  met  for  the  thermal 


instability,  the  Farley-Buneman  instability  can  be  strongly  ex¬ 
cited  in  much  shorter  time  scale  (millisecond).  This  is  because 
the  threshold  of  the  Farley-Buneman  instability  is  much  lower 
than  that  of  the  proposed  thermal  instability.  It  has  been 
shown  [Ossakotv  et  ai,  1975]  that  the  effects  of  finite  ku 
should  also  be  considered  in  determining  the  maximum 
growth  rate  of  the  Farley-Buneman  instability  for  such  large 
transverse  currents.  However,  we  don’t  expect,  in  reality,  that 
the  relative  drift  speed  can  be  too  much  larger  than  the  thresh¬ 
old  drift  of  the  thermal  instability.  Therefore  the  proposed 
instability  will  prefer  to  operate  for  the  case  of  minimum 
threshold  (fc,,  =0)  instead,  so  that  the  diffusion  damping  is 
minimized.  Though  the  proposed  instability  for  generating 
large  scale  irregularities  has  much  smaller  growth  rate  than 
that  of  the  Farley-Buneman  instability,  it  does  not  imply  that 
the  proposed  instability  will  not  have  significant  impact  re¬ 
garding  the  ionospheric  disturbances.  On  the  contrary,  the 
perturbations  introduced  by  the  large-scale  field-aligned  den¬ 
sity  irregularities  to  the  ionosphere  are  believed  to  be  impor¬ 
tant  by  the  facts  that  the  large  scale  field-aligned  density  ir¬ 
regularities  usually  have  high  saturation  levels  (e.g.,  fluctu¬ 
ations  in  a  few  percent  of  the  background  density  have  often 
been  observed  in  the  HF  heating  experiments),  and  can  stay 
for  a  very  long  period  of  time  even  after  the  driving  force  is 
removed.  While  our  proposed  thermal  instability  can  be  devel¬ 
oped  in  less  than  a  second,  E  x  B  instability  that  can  also 
excite  large  scale  irregularities  takes  a  few  tens  of  seconds  for 
its  development. 

The  excited  hundred  meter  scale  ionospheric  irregularities 
should  then  be  able  to  cause  spread  £  echoes  in  ionograms. 
Also,  in  the  absence  of  F  -egion  irregularities,  these  £  region 
irregularities  may  still  cau^e  amplitude  and  phase  scintillations 
of  beacon  satellite  signals.  These  phenomena  can  be  experi¬ 
mentally  studied  with  coordinated  incoherent  scatter  radar, 
ionosonde,  and  scintillation  measurement  techniques.  In  addi¬ 
tion  to  the  distinct  difference  in  scale  lengths,  the  irregularities 
induced  by  the  proposed  mechanism  are  purely  growing 
modes  (i.e..  zero  real  frequency),  while  those  excited  by  the 
Farley-Buneman  instability  have  real  frequencies. 
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Abstract 


An  instability  that  leads  to  the  filamentation  of  large-amplitude  Alfven  waves 
|  and  gives  rise  to  purely  growing  density  and  magnetic  field  fluctuations  is  studied. 

The  dispersion  relation  of  the  instability  is  derived,  from  which  the  threshold 
conditions  and  the  growth  rates  of  the  instability  are  analyzed  quantitatively  for 
|  applications  to  the  solar  wind  plasma.  We  have  examined  their  dependence  on  the 

filamentation  spectrum,  the  plasma  ft  and  the  pump  frequency  and  intensity  for  both 
right-hand  and  left-hand  circularly  polarized  Alfven  waves.  The  excitation  of 
filamentation  instability  for  certain  cases  of  interest  is  discussed  and  compared  with 
that  of  the  parametric  decay  and  modulation  instability.  The  relevance  of  the 
proposed  instability  with  some  observations  is  commented. 
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I,  Introduction 

A  great  deal  of  interest  in  the  stability  of  the  MHD  system  in  the  presence  of 
finite  amplitude  Alfv£n  waves  has  arisen  in  recent  years  [Wong  and  Goldstein, 
1986;  Terasawa  et  al.  1986;  Longtin  and  Sonnerup,  1986].  It  is  mainly  stimulated  by 
the  frequent  observations  of  large-amplitude  hydromagnetic  fluctuations  in  the  solar 
wind  at  1  AU  [Belcher  and  Davis,  1971],  in  the  high  speed  streams  of  the  solar  wind 
[Abraham-Shrauner  and  Feldman,  1977],  in  the  upstream  of  the  Jovian  bow  shock 
[Goldstein,  et  al.,  1985],  and  in  the  interplanetary  shocks  and  the  terrestrial  foreshock 
(Vinas,  et  al.  1984;  Smith,  et  al. ,  1985].  These  fluctuations  are  generally  believed  to 
be  circularly  polarized  Alfvfen  waves  propagating  almost  exactly  field  aligned. 
A  class  of  parametric  instabilities  excited  by  the  circularly  polarized  Alfv6n 
pumps  has  thus  been  investigated.  Two  different  instability  processes  are  currently 
discussed  in  the  literature.  One  leads  to  the  parametric  decay  instability  [Galeev  and 
Oraevskii,  1963;  Sagdeev  and  Galeev,  19669;  Cohen  and  Dewar,  1974;  Terasawa  et  al., 
1986],  and  the  other  one  gives  rise  to  the  modulation  instability  [Lashmore-Davies, 
1976;  Ionson  and  Ong,  1976;  Goldstein,  1978,  Derby,  1978;  Longtin  and  Sonnerup,  1986]. 
Moreover,  the  modulation  instability  has  also  been  analyzed  in  the  region  that 
describes  the  nonlinear  evolution  of  the  Alfv6n  waves  propagating  along  a 
static  magnetic  field  [Mio  et  al.  1976a  and  b;  Mjtflhus,  1976;  Spangler  and  Shering, 
1982],  A  derivative  non-linear  Schrodinger  equation  possessing  soliton  solutions  is 
derived  to  govern  the  evolution  of  nonlinear  Alfven  waves.  On  the  contrary,  Ovenden 
et  al.  [1983]  have  shown  that  the  evolution  of  nonlinear  Alfven  waves  is  governed  by 
a  se*  of  three  coupled  equations  which  in  turn  are  related  to  the  nonlinear 
Schrodinger  equation  with  known  soliton  solutions.  The  possible  applications  of 
Alfv6n  solitons  to  solar  and  astrophys ical  plasmas  have  been  discussed  by 
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Ovenden  et  al.  [1983]  and  Spangler  and  Sheerin  [1983]. 

Recently,  Sakai  and  Sonnerup  [1983]  investigated  the  effects  of  dispersion  on 
the  modulation  instability  excited  by  the  circularly  polarized  Alfven  waves.  Their 
analysis  is  restricted  to  the  situations  that  dispersive  effects  are  weak  and  the  wave 
number  k  of  the  sound  wave  is  much  smaller  than  the  wave  number  kQ  of  the  pump 
wave  (k  «  k^.  These  restrictions  have  then  been  lifted  in  the  more  recent  work 
by  Wong  and  Goldstein  [1986].  They  study  the  dispersive  effects  on  both  the 
modulation  and  decay  instabilities  in  a  unified  manner  over  a  wide  range  of  physical 
parameters.  In  addition,  the  results  show  that  the  maximum  growth  rate  of  the 
modulation  instability  occurs  as  k  is  comparable  to  kQ.  A  new  but  weaker 
instability  existing  in  a  very  narrow  bandwidth  near  k  «  kQ  has  also  been 
revealed  by  their  analysis. 

Ir.  this  paper,  a  filamentation  instability  which  has  not  been  considered  by 
the  previous  workers  is  analyzed.  It  is  known  that  a  large-amplitude,  initially 
uniform,  wave  propagating  in  a  plasma  can  break  up  into  filaments  because  of  the 
filamentation  instability  [G.  Schmidt,  1979;  Kuo  and  Schmidt;  1983].  This  starts  from 
small  perturbations  in  the  plasma  density  and  it  results  in  a  modulation  of  the  plasma 
dielectric  constant  and  wave  distribution,  which  in  turn  increases  the  density 
perturbations.  The  purpose  of  the  present  paper  is  to  show  that  the  filamentation 
instability  can  be  excited  together  with  the  parametric  decay  and  modulation 
instabilities  by  the  large-amplitude  circularly  polarized  Alfven  waves  in  the  solar  wind. 
The  threshold  conditions  and  the  growth  rates  of  the  instability  will  be  determined 
and  compared  with  those  of  the  decay  instability  and  the  modulation  instability.  Some 
observations  will  be  commented  for  corroborating  the  predicted  characteristic  of  the 
proposed  instability. 
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The  organization  of  the  paper  is  as  follows.  In  Section  II,  we  derive  the 
coupling  equations  for  Alfven  sidebands  and  the  purely  growing  magnetostatic  mode. 
A  dispersion  relation  is  obtained  in  Section  III  and  analyzed  for  the  threshold  fields 
and  the  growth  rates  for  various  cases.  The  numerical  results  are  also  given  in 
Section  III.  Finally  presented  in  Section  IV  are  a  summary  and  brief  discussions. 


4 


II.  Coupling  Equations 


We  investigate  the  propagation  and  fi lamentation  of  ducted  large-amplitude, 
circularly  polarized  Alfven  waves  in  infinite,  spatially  uniform  plasmas  embedded  in  a 

constant  magnetic  field  B0  *  zB0.  The  wave  magnetic  fields  are 
represented  by 


„  „  i(k  z-w  t) 

(x  ±  iy)B  e  °  0 

P 


+  c  .c . 


(1) 


where  Bp  is  the  unperturbed  wave  field  intensity,  assumed  to  be  constant  and  a 
real  quantity  for  simplicity,  the  ±  signs  refer  to  the  right-and  left-hand  circular 
polarizations,  respectively;  1^  and  wa  are  the  wave  number  and  angular  wave 

frequency,  satisfying  the  dispersion  relation  w02  *  (1  ± 

w0/q)k0IvAJ;  wherein  vA  and  are  the  AlfvSn  speed  and  the 

ion  cyclotron  frequency,  respectively. 

The  zeroth -order  velocity  responses  of  electrons  and  ions  to  the  Alfv£n  waves 
can  be  written  as 


v 

pet 


(1  ±lWVpi±  <* 


„  i(k  z-u)  t) 

iy) (w  Ik  )(B  /B  )e  °  ° 

o  o  p  o 


+  c.c. 


(2) 


where  the  subscripts  e  and  i  refer  to  electrons  and  ions. 


The  process  under  consideration  is  the  scattering  of  the  unperturbed  Alfven 
waves  into  sidebands  (propagating  oblique  to  the  z  direction)  by  the  simultaneously 
excited  density  perturbations  associated  with  the  purely  growing  magnetostatic  modes. 

Let  k  =  xk  be  the  filamentation  wave  vector,  the  density  and  magnetic 
field  perturbations  of  the  purely  growing  magentostatic  modes  have  the  expression 

ng  »  n^e^coskx 
and 

Bg  «  zBje^sinkx, 

where  n,  and  B,  are  real  amplitudes,  and  i  is  the  growth  rate; 
quasi  neutrality  has  been  assumed.  The  basic  equations  that  are  linearized  for 
analyzing  the  purely  growing  modes  include  the  continuity  equations,  the  momentum 
equations  for  both  electrons  and  ions,  and  the  Maxwell  equations: 

(3) 

(4) 

—  E  =  -  (y/c)B  and  —  B  =  (4im  e/c)(v  -v  )  (5) 

3x  sy  s  v  3x  s  o  sey  siy 

where  n0,  Te^  and  m(M)  are  the  unperturbed  plasma  density,  and  the 

unperturbed  electron  (ion)  temperature  and  mass,  respectively;  Fei  are 
the  nonlinear  Lorentz  forces  experienced  by  electrons  and  ions  that  reduce  to  the 


3  3 

vn  +  n  —  v  ~  0  ~  Yn  -hn  —  v 
1  s  o  3x  sex  s  o  3x  six 


-,mv  +  f  =  -  x(T  /n  )  ~  n  -  e(E  +  i  v  x  z  B  ) 
se  e  e  o  3x  s  s  c  se  o 


yMv  .  +  F.  =  -x(T./n  )  ^  u  +  e(E  +  -  v  xzB  ) 
si  i  i  o  3x  s  s  c  si  o 
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pondermotive  forces  in  the  unmagnetized  plasma  case. 


In  terms  of  v 


and  Svei±  that  represent  the  velocity  responses  of  electrons  and  ions  to 
the  unperturbed  Alfven  waves  and  the  sidebands  respectively,  the  nonlinear  Lorentz 
forces  are  given  by 


(6) 

where  the  notations  mti  -  m(M)  and  etj  =  +  1  are  used.  One  can 

show  in  a  self-consistent  way  that  the  nonlinear  Lorentz  forces  act  only  in  the  x 

direction,  i.e.,  Fg.  «  x  Fe-. 

We  now  substitute  the  expressions  vgg.  =  (x  v§e-x  +  y  v,ejy)e*sinkx,  vgex  = 

v«K  =  viix,  Eg  =y  ^^e^coskx  and  Fe-  =  x  Pe-e"rtsinkx  into  (3)-(5).  The  relation  Bg/Bc  = 

+  l^c2) ]¥f’g/n0  is  obtained  and  the  resultant  equations  can  be  combined  into  a 
single  coupled  mode  equation  for  the  purely  growing  magnetostatic  mode 


{-ik  6v  .  v 


o6ve,i±zVpe,i±X2"  Vpe.i±x[Vx(6Ve,i±X2)1}  +  c‘c' 


7 


{T2  +  [I «*vA*/(  1  +  k*v»A/f^H))  +  k2c2s)n./n0  -  (k/M)(F,  ♦  F,) 


(7) 


This  equation  gives  the  purely  growing  density  perturbation  as  a  function  of  the 
Alfven  waves  and  sidebands.  Without  the  nonlinear  Lorentz  forces,  Eq.  (6)  reduces  to 
the  linear  dispersion  relation  of  the  magnetosonic  eigen-mode.  The  magnetostatic 
mode  is,  however,  a  nonli nearly  driven  mode. 

The  Alfven  sidebands  are  excited  through  the  beating  current  density  driven  by 
the  Alfven  pump  wave  fields  on  the  density  perturbation  of  the  magnetostatic  mode. 
The  coupled  mode  equation  for  these  sidebands  can  be  derived  from  the  following 
fluid  equations  for  electrons  and  ions  together  with  the  Maxwell  equations: 


<5n  +  V*(n  6v  +  n  v  )  =  0  *  — -  6n  +V-(n  6v  +  n  v  ..) 

3t  ±  o  e±  s  pe±  3t  t  o  i±  s  pi± 


6E 


+  —  6v  x  zB  *  -  (T  /n  e)  V  6n 
±  c  e±  o  e  o  ± 


(8) 

(9) 


6v .  =-c^V(6n  In  )  +  (e/Mc)  (6v  .  -  5v  )  *  zB  (9  ’) 

3t  i±  s  ±  o  i±  e±  o 

Vxil  =  --  (10) 

±  C  dt  ± 


and 


Vx6^±  =  (4ivnoe/c)[(6vi± -6ve±)  +  (iWn^  ^pi± _V± 


)1 


(ID 


where  c  = 
s 


(Te  +  Ti)/mi- 


g 


We  have  neglected  the  nonlinear  Lorentz  forces  ("v  »Vv)  in  (9) 
and  kept  the  beating  currents  (ntvp+)  in  (8)  and  (11)  as  the  driving 
sources  of  the  Alfven  sidebands.  The  underlying  reason  is  that  the  purely  growing 
mode  is  more  effective  in  producing  density  perturbations  than  velocity  perturbations. 

Following  the  functional  dependence  of  the  magnetostatic  mode,  we  express  the 
physical  quantities  of  the  sidebands  as 


6n± 


,  -  -  -  ~  -  ►  l(k  Z-  W 

[  (x  6v  ±  iy  6v  )coskx  +  i  z  <5v  sinkxle^  e  °  °  +  c.c. 

eitx  ei±y  ei±z 


i(k  z-o)  t) 

i  6n+sinkx  e'  e  +  c.c. 


i(k  z-w  t) 

[ (x  t ia+y)coskx-lz(k/ko)sinkx]6B+  e1  e  °  0  +  c.c.  (12) 


where  the  upper  sideband  (k+  =  kQ  +  k,  w+  =  w0  +  n)  and  lower  sideband  (k.  =  kQ-k, 
w.  =  w0  -(n)*)  propagate  together  along  the  magnetic  field  and  form  a 
standing  wave  pattern  across  the  magnetic  field,  a+  are  unknown  parameters 
defining  the  polarizations  of  the  Alfven  sidebands  and  will  be  determined  shortly;  the 
expressions  of  6E+  can  be  obtained  from  the  Faraday  law  (10). 

Substituting  (12)  into  ( 8) -( 1 1) ,  and  eliminating  5n+  from  (8)  and  (9), 


we  first  obtain 


5v  =  -  (lj  /k  )5B  /B 
etx  o  o  to 


5v 


eiy 


(u>  /k  )6B  /B 


ZOO 


6v .  =  t(k  vj/u  )(n  /n  )  (B  /B  )-(k  ftj/k2w  )(u2/kV4-l)(l?a  u  /fl.) 

i±x  oAo  so  po  oi  o  ooti  ±oi 


(6B  /B  ) ] / f 
±  o 


«»Ily  *  *<v„>  l(vX><VV(VV  +  (“o/to)[l  +  (Vk>2<V“o) 

(“o/kovti  •  1)<*ai"“o/r,l)  1  (iEr/Bo)l/f 

5vUz  ~<ko/k>  (^0>  <%'-„>  <  V»„)  +  <“o/ko> 

(i+  a+fW^)  (6B+/Bo)  }/f  (13) 


where  f  =  1  +  (k0/k)2(l\/wo)2(l-w0J/^2)(w02/k02vu-l),  and  vti  =  (Tj/M)1/2  is  the  ion 
thermal  speed. 

The  missing  component  5ve+I  can  be  obtained  from  the  z 
component  of  (11),  it  reads 


6v 

etz 


(14) 
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These  relations  of  (13)  and  (14)  are  then  used  in  the  remaining  two  equations, 
i.e.,  the  x  and  y  component  of  (11),  leading  to  the  determination  of  the  parameter 

a+  =  {(l±o)o/fii)ll-(ko/k)2(a)J/kVi-l)]  +  (f  ifl^Kk^/flJxi  +  k2/^)}/ 

{ (1/f  ±(Jo/ni)[l-(ko/k)2(^/k2vt2-l)l 
+  (i-n2/u2f)[fk2v2/n2-(ko/k)2(co2/k2v2ri)J}  (15) 

and  that  of  the  coupled  mode  equations  for  the  Alfvdn  sidebands 

U-(ko/k)2(n1/Uo)2(i+«±»0/al)(^/k^1-i)/f  i  »i(k2v2/nt»o)K«Bt/»0) 

-  -U/f  iUo/ai)(k2v2/^)(is/„o)(Bp/Bo)  (16) 

Equation  (16)  shows  the  sideband  fields  to  be  produced  by  the  pump  field  and 
the  purely  growing  density  perturbation  of  the  magnetostatic  mode.  Thus,  equations 
(7)  and  (16)  form  a  complete  coupled  set  of  equations  describing  the  proposed 


instability  process. 


III.  Dispersion  Relation 


With  the  aid  of  (2),  (13),  (14)  and  (6),  the  coupling  term  on  the  right  hand  side 


of  (7)  can  be  expressed  explicitly  as 


(Fe  +  F.)/M  =  ±2k(n./a,o)(kov^i)2({i+k2/k2  _  (a,oQ1/k2v2)[a±+(n/M) 


(i  ) (i  ± w  /n  ) 

±  o  i 

(1/f )  (ko/k)2(r_/too)2(l  +  Ct±  “0/ni))}(6B±/Bo)  +  (Q^w  ) 


{(i  ±2w  /a.) 
o  i 


(17) 


Substituting  (17)  into  (7)  and  combining  the  resultant  with  (16),  we  finally  obtain 
the  dispersion  relation 


Y2  +  k2v2/(l  +  k2v2/n  G.)  +  k2c2  +  2k2(k  v2/u>  )2[(1  ±  2u)  /G.)-(k  /k)2 
A  Aei  s  oAo  olo 

*  -2k2(V*/“o)2(kov*/Bi)(1  =  V“of)U+k2/ko  +  <“„VkovA> 

[a++  (tn/M)  (1  +  a  +  )  (1 

-(l/f)(ko/k)2(Gi/ujo)2(l  +  a+Wo/ni)]}(Bp/Bo)2/ 

u-(ko/k)2(Ul/-o)2(i  *«!"0/v<“2/ko''t2rl>/E  a8) 
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Equation  ( 18)  is  a  general  expression  of  the  dispersion  relation  for  fi lamentation 
instability  of  ducted  Alfven  waves.  It  will  be  analyzed  in  the  following  to  determine 
the  threshold  conditions  and  the  growth  rates  of  the  instability. 

In  principal,  (18)  is  a  quadratic  equation  in  n  and  can  be  solved  analytically. 
However,  since  the  coefficients  of  (18)  are  complicated  functions  of  the  constant 
parameters  0  **  c*,/v*A  and  K.  =  K0VA/fy  and  the  normalized 
variable  parameter  K  ~  k/k0  notations  used  in  Wong  and  Goldstein  [1986]  are 
adopted).  A  numerical  analysis  of  (18)  with  different  parameters  is  desirable.  In 
doing  so,  we  vary  K.  from  0  to  2,  and  examine  the  dependences  of  the  threshold 
intensity  and  the  growth  rate  on  0  and  K  for  both  right-(R)  and  left-hand  (L-H) 
circularly  polarized  Alfven  waves.  Presented  in  Figures  1-4  are  the  dependence  of  the 
normalized  threshold  intensity  (f?*th)  on  K.  17 2th 

( Bp/ Bo)2thi  s  obtained  by  setting  *7  =  0  in  (18).  Figures  1  and  2 

correspond  to  the  R-H  circularly  polarized  waves.  It  is  shown  that  the  filamentation 
instability  has  a  lower  threshold  level  for  a  wave  with  a  smaller  frequency  (i.e. 
smaller  ic).  In  addition,  the  threshold  of  the  instability  also  varies  with  0.  For 
very  low  frequency  waves  (*  <  0.05)  the  threshold  intensity  increases  with  0. 
It,  however,  becomes  a  decreasing  function  of  0  for  it  >  0.1.  This  trend  has 
been  demonstrated  in  these  two  Figures,  wherein  k  =  0.01  and  0.3  are  used  as  the 
representative  parametric  values.  We  have  excluded  the  threshold  values  in  the 
neighborhood  of  X  (marked  on  the  curves)  from  the  Figures.  This  is  because,  in  that 
wave  number  region  the  nonlinear  Lorentz  forces  on  electrons  and  on  ions  are  in 
opposite  direction  with  comparable  magnitude  so  that  the  two  terms  proportional  to 

(Bp/B0)2  in  (18)  tend  to  cancel  to  each  other  and  the  threshold,  thus,  becomes 
very  high. 
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Figures  3  and  4  correspond  to  L-H  circularly  polarized  waves.  It  is  found  that 
the  instability  can  only  be  excited  by  waves  with  higher  frequencies  (e.g.  *  > 
0.15)  with  reasonable  thresholds.  The  dependence  of  the  threshold  on  k  for  p  =  1.5 
and  that  on  P  for  «  =  0. 3  are  also  shown  in  Figures  3  and  4,  respectively.  It 
appears  that  the  instability  prefers  to  be  excited  in  the  region  K  <0.5.  Otherwise, 
the  threshold  increases  very  fast  with  K.  In  the  region  K  <  0.5,  the  threshold 
decreases  with  both  /c  and  P  until  the  instability  becomes  forbidden  for  the  region 
of  very  small  <c  (e.g.  the  region  below  K.  =  X  marked  in  both  Figures  3  and  4  for 
k  =  0.3  and  P  =  1.5). 


We  next  present  the  growth  rates  of  the  instability.  Shown  in  Figures  5  and  6 
are  the  functional  dependence  of  the  growth  rate  7  on  the  wave  field  intensity  17* 
of  R-H  wave  for  low  0.5)  and  high  fl(=  1.5)  cases,  where  k  and  K  are 
considered  to  be  constant  parameters.  In  order  to  obtain  adequate  information  for 
the  dependencies  of  7  in  these  two  cases,  two  representative  values  for  each  of  the 
two  parameters  *  and  K  are  chosen  in  the  figure  presentation,  /c  =  0.01  and  0.3 
stand  for  the  cases  of  the  low-frequency  and  high-frequency  waves,  respectively.  We 
then  use  K  =  0. 1  and  1  to  characterize  the  regions  of  large-scale  and  small-scale 
filamentation  instability,  respectively.  The  results  show  that,  in  general,  the  growth 
rate  increases  with  P  and  K  and  it  decreases  with  k.  For  L-H  wave  case,  the 
instability  can  only  be  excited  by  the  high-frequency  waves  in  the  high  P  plasma. 
We,  therefore,  choose  P  =  1.5  and  *  =  0. 15  and  0.3  to  evaluate  the  dependence  of 
7  on  rf.  Again  K  =  0. 1  and  1  are  considered  and  the  results  are  presented  in 
Figure  7.  In  this  case,  the  growth  rate  increases  with  «  for  K  =  0. 1  and  becomes  a 
decreasing  function  of  *  for  K.  =  1. 
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IV.  Summary  and  Discussion 

We  have  investigated  the  filamentation  instabilities  of  large  amplitude,  circularly 
polarized  Alfven  waves  propagating  along  the  background  magnetic  field.  The 
instabilities  are  excited  via  the  scattering  of  the  unperturbed  Alfven  pumps 
(considered  as  a  pump  wave)  into  sidebands  by  the  density  perturbations  that  are 
associated  with  the  simultaneously  excited  purely  growing  magnetostatic  modes.  A 
four- wave  coupling  process  is  then  considered  for  the  analysis  of  the  instabilities. 
The  theory  developed  is  based  on  the  two  fluid  plasma  model.  In  general,  the  fluid 

model  is  valid  when  the  condition  k^/nf  «  K2k?£/2  «  1  is  satisfied.  In  the 

present  work,  the  maximum  value  of  the  parameter  K2k?0/2  for  K  =2,  k  =  0.3 
and  P  =  1.5  is  0.27,  which  reasonably  justifies  the  use  of  the  fluid  model.  For  the 
more  general  consideration  that  includes  the  region  K2k0/2  >  1,  however, 
a  kinetic  plasma  model  should  be  used. 

The  nonlinear  source  for  the  Alfven  sidebands  is  the  beating  current  driven  by 
the  pump  wave  field  on  the  density  perturbation  of  the  purely  growing  mode,  whereas, 
the  nonlinear  Lorentz  force  (which  reduces  to  the  pondermotive  force  in  the 
unmagnetized  case)  introduced  by  the  spatial  gradient  of  the  resultant  high  frequency 
wave  field  is  the  driving  source  for  the  nonoscillatory  (i.e.,  the  purely  growing)  mode. 
These  nonlinear  effects  result  in  the  coupling  of  nonoscillatory  mode  with  the  Alfven 
sidebands  through  the  Alfven  pump  wave.  A  dispersion  relation  of  the  instabilities  is 
thus  derived  by  combining  the  coupled  mode  equations  together.  Solving  the 
dispersion  relation,  we  have  determined  the  threshold  fields  and  the  growth  rates  for 
the  caes  of  the  R-H  and  L-H  circularly  polarized  pump  waves.  We  have  also 
determined  the  functional  dependencies  of  the  threshold  field  »?th  on  the 
wavenumber  K  of  the  nonoscillatory  mode,  the  0  of  the  plasma,  and  the  frequency  * 
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of  the  pump  (see  Figures  1-4).  The  dependence  of  the  growth  rate  7  on  *he  pump 
intensity  r?  is  examined  for  several  sets  of  representative  parameters  (K,ftic) 
(see  Figures  5-7). 

For  the  geophysical  applications  of  the  concerned  instabilities,  we  follow  the 
examples  discussed  by  Wong  and  Goldstein  [1986],  The  finite  amplitude  Alfven  waves 
that  were  observed  in  the  high  speed  streams  of  the  solar  wind  and  discussed  by 
Abraham-Shrauner  and  Feldman  [1977]  are  the  left-hand  circularly  polarized  with 
k  *  0. 3  and  t?  =  0.05.  As  shown  in  Wong  and  Goldstein  [1986],  these  waves  are 
stable  to  both  the  modulational  and  the  decay  instabilities  for  P  >  1.  From  our 
analysis,  however,  it  is  found  that  these  waves  are  unstable  to  the  fi lamentation 
instability  in  the  small  K.  region.  The  growth  rates  of  the  instability  evaluated  for 
P  =  1.5  have  been  presented  in  Figure  8.  We,  therefore,  predict  the  appearance  of 
nonoscillatory,  cross-field  plasma  density  perturbations  and  magnetic  field 
perturbations  in  the  background  solar  wind  plasma.  However,  these  predicted  cross- 
field  density  striations  and  magnetostatic  structure  probably  have  not  been  observed 
in  Abraham-Shrauner  and  Feldman  [1977],  The  difficulty  in  the  in-situ  measurements 
by  satellites  lies  on  the  large  velocity  difference  between  the  solar  wind  and  the 
satellite. 

Low  frequency,  left-hand  polarized  finite  amplitude  Alfven  waves  have  also  been 
observed  in  the  upstream  of  the  Jovian  bow  shock  [Goldstein,  et  al.,  1985],  It  is 
shown  from  our  analysis  that  this  wave  is  filamentally  stable.  In  fact,  only  the  decay 
instability  is  very  weakly  unstable  for  the  case  k  =  0.02,  rp  =  0.05  and  P  =  1.5. 
Whereas,  the  observation  of  right-hand  polarized  waves  were  not  reported.  It  can  be 
speculated  from  the  linear  dispersion  relations  of  the  Alfven  waves.  Because,  in  such 
low  frequency  region,  the  dispersion  relations  of  right  and  left-hand  polarized  Alfven 
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Figure  Captions 


Fig.  i 

Fig.  2 

Fig.  3 

Fig.  4 

Fig.  5 

Fig.  6 

Fig.  7 

Fig.  8 

Fig.  9 


Functional  dependence  of  the  threshold  field  intensity  on  K  and  P  for  low 
frequency  (<c  =  0.01)  R-H  circularly  polarized  pump. 

Dependence  of  the  threshold  field  intensity  on  K  and  P  for  high  frequency 
(k  =  0.3)  R-H  circularly  polarized  pump. 

Dependence  of  the  threshold  field  intensity  on  K  and  k  for  L-H  circularly 
polarized  pump  in  high  #=  1.5)  plasma. 

Dependence  of  the  threshold  field  intensity  on  K  and  P  for  high  frequency 
(ic  *=  0.3)  L-H  circularly  polarized  pump. 

The  dependence  of  the  normalized  growth  rate  7  =  7/k„VA  on  the 
R-H  circularly  polarized  pump  intensity  if,  where  K  =  0. 1  and  (*,  P) 
as  indicated  are  set  in  the  evaluation. 

Dependence  of  7  on  if  of  R-H  circularly  polarized  pump  for  K.  =  1 
and  (k,  P)  as  indicated. 

Dependence  of  7  on  rf  of  a  L-H  circularly  polarized  pump  in  a 
plasma  with  P  =  1.5. 

Dependence  of  7  on  K.  for  a  L-H  circularly  polarized  pump  with 
=  0.05  and  k  =  0.3  in  a  plasma  with  P  =  1.5. 

Dependence  of  7  of  K  for  a  R-H  circularly  polarized  pump  with  if 
=  0.05  and  k  =  0.02  in  a  plasma  with  P  =  1.5. 


21 


Functional  Dependencies  of  the  threshold  field  intensity  on  K  and 
/5 for  low  frequency  (*  =  0.01)  R-H  circularly  polarized  pump 
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dependence  of  the  normalized  growth  rate  /= //k0VA  on  the 
circularly  polorized  pump  intensity^2,  where  K=0.1  and 
as  indicated  are  set  in  the  evaluation. 
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Dependence  of  Y  on  T^of  a  L-H  circularly  polorized  pump 
in  a  plasma  with  1.5 


2.23 


4.5 


Dependence  of  / of  K  for  a  R-H  circularly  polarized  pump 
w jth  T)2  =0.05  and  *  =  0.02  in  a  plasma  with  /?  =  1.5. 
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Abstract 


This  work  is  intended  to  establish  a  useful  method  for  studying  the  coupling  and 
propagation  of  hydromagnetic  waves  in  the  dipole  model  of  the  magnetosphere.  A 
numerical  scheme  for  solving  a  system  of  two  mixed-type  coupled  partial  differential 
equations  having  variable  coefficients  is  developed.  The  stability  condition  of  the 
scheme  has  been  examined  by  the  Fourier  method  and  matrix  method.  An  initial  value 
problem  has  then  been  solved  to  demonstrate  the  relevance  and  applicability  of  the 
work  to  the  study  of  magnetospheric  hydromagnetic  waves. 
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Introduction 


In  the  study  of  wave  propagation  in  inhomogeneous  magnetoplasma,  the  physical 
system  is  generally  characterized  by  coupled  partial  differential  equations.  One 
specific  example  is  the  geomagnetic  micropulsation  originated  by  the  hydromagnetic 
waves  in  the  magnetosphere  [1-3].  It  is  shown  that  these  waves  in  the  poloidal  and 
toroidal  modes  are  generally  governed  by  two  coupled  second  order  partial  differential 
equations  [4,  5],  Though  these  two  equations  can  be  reduced  and  combined  into  a 
single  second  order  ordinary  differential  equation  for  a  simplified  cylindrical  model  of 
the  magnetosphere  [5],  there  is  no  apparent  way  to  simplify  the  system  of  coupled 
equations  for  a  more  realistic  dipole  model  of  the  magnetosphere  [4  ]  in 
which  a  dipole  geomagnetic  field  is  assumed.  Moreover,  the  equation  for  the  toroidal 
mode  is  shown  to  be  a  parabolic  type  of  partial  differential  equation,  while  the 
poloidal  mode  is  elliptic.  Both  equations  have  variable  coefficients,  therefore,  for 
practical  applications,  there  is  a  considerable  interest  in  finding  or  developing  means 
of  dealing  with  such  a  system  of  two  mixed-type  coupled  partial  differential  equations 
with  variable  coefficients. 

Because  of  the  difficulty  anticipated  and  approximations  necessary  in  the 
analytical  methods  [6,  7],  a  numerical  approach  is,  therefore,  thought  to  be  more 
feasibile.  Various  numerical  methods  have  been  developed  for  dealing  with  partial 
differential  equations.  However,  the  applicability  of  each  method  is  mainly  determined 
by  its  range  of  stability.  As  early  as  the  1960’s  the  energy  method  had  been  used  to 
solve  certain  classes  of  the  partial  differential  equations  with  variable  coefficients  [8, 
9].  The  range  of  applicability  of  the  method  has  also  been  discussed  and 
demonstrated  for  limited  cases.  Thus,  the  elliptic  type  partial  differential  equation 
with  varible  coefficients  has  to  be  solved  by  the  high-order  methods  [10,  II],  On  the 
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other  hand,  the  Galerkin  method  [12]  has  been  shown  to  be  more  applicable  for 
dealing  with  an  equation  of  the  form  ut  =  P(x,  D)u,  where  D  is  a  differential 
operator  with  respect  to  x.  It  is  noted  that  all  the  above  mentioned  methods  have 
been  developed  based  on  the  finite  difference  techniques. 

The  Fourier  method  is  widely  used  for  the  stability  analyses  of  the  finite 
difference  schemes  [13,  14].  However,  an  analytical  study  of  the  stability  of  the 
difference  schemes  by  the  Fourier  method  was  often  found  to  be  very  complicated  or 
practically  impossible.  Numerical  approaches  were  thus  introduced  [15,  16].  In 
addition  to  the  Fourier  method,  a  matrix  method  [17,  18]  which  includes  the  boundary 
conditions  in  the  analysis  was  also  suggested  for  the  stability  analysis.  In  the  present 
work,  both  method  will  be  employed  in  the  stability  analyses.  The  Fourier  method 
will  first  be  used  to  examine  parametrically  the  stability  conditions  of  the  designed 
scheme.  These  conditions  will  then  be  incorporated  in  applying  the  matrix  method 
which  is  thus  used  to  exemplify  the  validity  of  these  stability  conditions. 

For  the  Fourier  method,  solutions  of  the  monochromatic  wave  form  are 
substituted  into  the  difference  scheme  and  the  stability  boundaries  are  obtained  by 
von  Neumann  condition.  In  the  second  approach,  the  difference  equations  will  be 
arranged  into  a  matrix  form  AVj+1  =  BVj,  where  A  and  B  are  matrices  and 
Vj+1  and  Vj  are  solution  vectors  at  grid  number  j  +  1  and  j.  The  stability 
of  the  finite  difference  approximation  is  then  ensured  if  the  absolute  values  of  all  the 

eigenvalues  of  A_IB  are  less  than  or  equal  to  1. 

In  Section  II,  the  coupled  hydromagnetic  wave  equations  in  the  dipole  coordinate 
system  are  introduced.  They  represent  a  set  of  mixed-type  coupled  partial  differential 
equations  with  variable  coefficients.  The  finite  difference  equations  which  are  derived 
by  a  designed  numerical  scheme  for  the  system  of  partial  differential  equations  are 
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presented  in  Section  III,  in  which  the  stability  analyses  by  Fourier  and  matrix  method 
are  also  presented.  As  a  specific  example,  the  boundary  conditions  of  the  wave 
equations  are  specified  in  Section  IV.  Presented  in  Section  V  are  the  stability 
conditions  of  the  numerical  scheme  and  the  numerical  solutions  of  the  wave  equations. 
Finally,  a  summary  and  conclusion  are  in  Section  VI. 


II.  Governing  Equations 


Considered  next  will  be  the  propagation  and  coupling  of  hydromagnetic  waves  in 
the  magnetosphere.  The  coupled  hydromagnetic  wave  equations  in  the  dipole  model  of 
the  magnetosphere  are  given  by  [4] 


HiH2  I?  ♦  HtG2  |  +  (^-m2H1)«l/  =  mH,  |  ^  (I) 

H1H2  df. ?  £ *  +  H2G1  +  ^  ^2  J  ~  "^2  (*  “mI^2  &  €v  ^ 

where  (v,  fi,  <f>)  form  an  orthogonal  set  of  the  dipole  coordinates;  and 
are  the  normalized  toroidal  (shear)  and  poloidal  (compressional)  mode  fields  of 
the  hydromagnetic  waves. 


ii  _  _]_  u  _  K*'  +  Vr3) 
1  _  it3  ’  2  ~  r 


5H,  3H- 

G1  =  TjT  ’  °2  =  -5T 


dH 


r  is  related  to  v  and  n  through  the  equation  w  +  /i2r4  =  1, 


A2  =  A2  ( l+3/izr4)  (r)n 


-6 


Bo 

where  An  =  -==  is  the  Alfven  speed  at  the  equator  of  the  earth,  a  mode!  plasma 
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mass  density  profile  p  =  pQ(r) the  plane  of  equator  is  assumed;  B0  is 
the  geomagentic  field,  Re  is  the  earth’s  radius,  m  is  the  azimuthal  (i.e.  east-west 
direction)  mode  number  and  w  is  eigenfrequency  of  waves. 

Equations  (1)  and  (2)  constitute  a  set  of  mixed-type  coupled  partial  differential 
equations  with  variable  coefficients.  Equation  (1)  is  a  parabolic  type  of  P.D.E.,  while 
Eq.  (2)  is  an  elliptic  type.  For  the  convenience  of  the  following  analysis.  Equations 
(1)  and  (2)  are  rewritten  as 


=  b1(/*,v)Pw4 
where  P 

ai 
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XM  =  d[l*'  **  =  duX'  X  =  P’T’  al’  a2’  3S 

There  is  no  apparent  way  to  solve  (3)  and  (4)  analytically.  A  numerical  approach  is 
therefore  employed  in  the  present  study. 

III.  Finite  difference  discretization  and  stability  analyses 

The  governing  Eq.  (3)  and  (4)  are  approximated  on  a  non-staggered  grid  with 
grid  increments  Am  and  A v.  Both  the  computational  grids  used  are  carefully 
chosen  to  ensure  an  algebraic  system  of  equations  that  reach  second  order  accuracy. 
The  standard  two  and  three-point  centered  difference  formulas  to  approximate  first 
and  second  derivatives  with  respect  to  m  and  v  are  applied  in  this  finite  difference 
discretization.  Moreover,  a  parameter  0  valued  in  the  range  (0,  1)  is  introduced  as  a 
weighting  factor  into  the  system  of  difference  equations.  Thus,  the  numerical  scheme 
of  the  difference  equations  can  change  from  fully  implicit  (0  =  1)  to  fully  explicit 
(0  =  0).  Based  on  the  above  discussion,  the  resulting  finite  difference  equations 
converted  from  (3)  and  (4)  become 
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where  i:  1  — »  I  - 1  and  j  :  1  — ►  J  - 1 

Two  methods  are  used  to  establish  the  stability  conditions  of  the  finite 
difference  scheme  which  converts  the  differential  equations  (3)  and  (4)  into  the  finite 
difference  equations  (5)  and  (6). 

A.  Fourier  method 

In  this  method,  the  independent  solutions  of  the  difference  equations  are  all 

of  the  forms  pj  =  £  e'k^A*1  and  Tj  =  e'k/Ap  where  k 
is  a  real  spatial  wave  number  which  can  have  any  value.  Using  these  expressions  for 
pj  and  T j  in  (5)  and  (6),  yields 

aJ+1  aJ+1  aJ  aJ 

P  +  B  FjT  =  P  +  ( 1-0) FjT  (7) 

aJ+1  A-i+1  AJ  aJ 

e  F2P  +  F3f  =  (l-0)F2P  +  F«f  (8) 
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where  F,  =  £  (coskAp-1)  +  2a3A*  +  i(2a2  ^  sinkAp)  “] 

F2  =  "  bi^2  (coskAp-l)  +  2b3Ai/  +  i(2b2^j  sinkA/x) 


F3  {( l-a3b4)-aIb4  — 3{coskA/i-I)-^a11/b4  ^^(coskA/i- 1)+  2a3l/b4A p] 


sinkA/i 


^2b<  +  2a2,b<  f^inkAp)} 


F4  =  {( l-a3b4)-ajb4  ,  (coskA/j-1)  +  (l-0)[a1(,b4  ^^(coskA^-])  +  a3l/b4Aj-'] 


1  ta2b«  S'  -2(  l-5)a2l<,b4^sinkA/i]) 


The  coefficients  3j,  a2,  a3,  a^,  a2l/,  a3l,,  bj,  b2  and  b3  in  the  above 
expressions  are  evaluated  at  /1  =  / An  and  v  =  jA^.  Eq.  (7)  and  (8)  can  then  be 


expressed  in  the  matrix  form  as 


n  ^1 

L  *F=  Fs  J 


I  0-*)Fj 

0-«)F,  F . 
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which  leads  the  amplification  matrix 


G  = 


-  Fg-fld-flFjFa 
Fs-^f, 


0-2J)F, 

_  Fj-^F, 


( l-^)F1F3-ffr1F4  ~ 
F3-fl2F1F2 

F4-fl[l-g)F1F2 
F3-e2F2F2  _ 


The  definition  of  stability  here  is  the  classical  von  Neumann  stability  condition 
which  requires  that  the  norm  of  the  amplification  matrix  have  to  be  less  than  or 
equal  to  1.  It  is  noted  that  the  amplification  matrix  has  to  be  determined  for  every 
grid  point  because  of  the  variable  coefficients  appeared  in  the  coupled  equations  of 
the  present  study.  The  norm  of  the  amplification  matrix  at  each  grid  point  is  then 
I  evaluated  numerically.  The  result  of  a  parametric  stability  analysis  obtained  by 

varying  the  parameters  Ap,  Ai/  and  9  will  be  presented  in  Section  II  evaluated 
numerically. 

I 

B.  Matrix  method: 

By  transforming  the  difference  equations  (5)  and  (6)  into  a  heptagiagonal 
matrix  form,  it  reads 

Ai+1Ej+1  =  BjlEjl  (9) 

where 
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Matrices  A  and  B  have  dimension  2(1-1)  x2(  I  - 1 ) ,  while  matrix  E  has  dimension  2(  I  - 1 )  x  1 . 


The  stability  of  the  scheme  is  then  governed  by  the  norm  of  the  matrix  A_1B 


V.  Boundary  conditions 


In  the  magnetosphere,  the  inner  boundary  is  located  at  the  outside  edge  of 


the  ionosphere  (r  =  2Re),  while  the  outside  boundary  is  at  the  magnetopause  (r  = 


10Re).  We  thus  have  the  inner  boundary  at  v  =  0.5  and  outer  boundary  at  v  - 


0.1.  The  other  two  boundaries  are  the  north  pole  (m  =  0.25)  and  south  pole  (m  = 


-0.25).  These  boundaries  define  the  region  -0.25  <  M  <  0. 25  and  0. 1  <  v  <  0. 5  for 


the  magentospherc. 
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Assuming  that  the  ionosphere  is  a  perfect  conductor,  hence,  the  tangential 
components  of  electric  field  at  the  inner  boundary  are  equal  to  zero,  which  lead  to 
the  fixed  boundary  conditions:  P(u,  n  =  +  0.25)  =  0  =  T(x/,  jx  =  ±  0.25).  The 
boundary  conditions  at  the  magnetopause  =  0.1,  n)  and  T (v  =  0.1,  /x)  are,  in 

general,  imposed  by  the  external  source  such  as  the  solar  wind.  They  can  thus  be 

assumed  to  be  any  funcitons  of  /x  for  different  physical  situation.  In  the  present 
work,  we  set  the  boundary  conditions  at  u  =  0. 1  to  be 
P(0.1,  n)  =  P e(/x)  =  |sin(4irh*x)| 

T(0. 1,  n)  =  0 

where  PQ(p>  is  a  spatially  distributed  sinusoidal  perturbation;  h  is  the  harmonic 
number  of  sinusoidal  perturbation;  T  =  0  is  for  the  practical  reason  since  the 
perturbations  in  the  frequency  range  of  the  toroidal  mode  at  the  outer  boundary 
cannot  propagate  into  the  magnetosphere. 

The  frequency  w  is  given  to  be  2*  x  0.002Ht  by  the  experimental  data  [19], 
Two  extreme  cases  with  the  azimuthal  wave  number  m  =  20  and  100  respectively  will 
be  considered  consistent  with  the  observations  (|m|  =  20  -100)  (19-22).  We 
will  also  assume  a  reasonable  density  index  n  =  3  in  the  numerical  analysis.  The 
stability  conditions  of  the  numerical  scheme  and  numerical  solutions  will  be 
determined  by  using  these  boundary  conditions  and  parameter  values. 

VI.  Results 

In  order  to  find  the  stability  conditions  by  Fourier  method  which  is 
described  in  Section  III  A,  the  eigenvalues  of  amplification  matrix  G  are  evaluated  at 
every  grid  point.  By  varying  the  parameters  A/x,  Ax'  and  9,  the  roots  of  a 
second  order  equation  defined  by  the  eigenvalue  equation  are  determined  iteratively  to 
satisfy  the  von  Newmann  stability  condition  which  requires  the  absolute  values  of  both 
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roots  to  be  less  than  or  equal  to  1.  Fig.  1  shows  the  y  dependence  of  the  maximum 
absolute  value  of  both  roots  of  the  second  order  equation  over  the  entire  range  [- 
0.25,  0.25]  of  M  for  several  different  8  values  and  for  the  case  m  =  100,  where 
A/i  =  0.004  and  Ay  =  0.01334  are  chosen  for  the  reason  stated  in  the  following. 
It  should  be  pointed  out  that  only  8  and  At/  are  the  sensitive  factors  in  the 
scheme.  Because  the  grid  number  in  the  /i  direction  will  be  the  dimension  of  the 
matrices  A  and  B  of  Eq.  (9)  used  for  the  analysis  of  the  matrix  method.  Hence, 
A/*,  the  grid  size  in  the  M  direction,  is  limited  by  the  memory  size  of  the 
computer,  and  determined  reasonably  to  be  0.004.  In  other  words,  the  dimension  of 
the  matrices  A  and  B  is  set  to  be  248  x  248  and  the  program  is  executed  in  the  Cray 
X-MP/48  of  the  Pittsburgh  supercomputing  center.  As  shown  in  Fig.  1,  the  larger  8 
is  used  in  the  scheme,  the  more  stable  scheme  can  be  found.  Therefore,  the  fully 
implicit  scheme  (8  =  1)  will  be  chosen  to  be  the  second  parameter.  For  the  accuracy 
of  the  numerical  solution,  there  is  a  need  to  find  the  minimum  A y  in  the  numerical 
scheme.  By  specifying  A/i  =  0.004  and  0=1,  the  functional  dependence  of  the 
maximum  eigenvalue  of  amplification  matrix  G  on  v  for  different  Ay  values  is 
shown  in  Fig.  2.  From  this  analysis,  the  minimum  Ay  =  0.01334  is  determined.  We 
next  apply  the  matrix  method  to  examine  the  stability  conditions  obtained  by  the 

Fourier  method.  As  shown  in  Fig.  3,  all  of  the  eigenvalues  of  the  matrix  A_1B 
in  Eq.  (9)  are  found  to  be  less  than  1.  Therefore,  the  designed  numerical  scheme  has 
been  proven  to  be  stable  consistently  by  Fourier  and  matrix  method.  In  other  words, 
a  viable  numerical  algorithm  dealing  with  the  mixed-type  coupled  partial  differential 
equations  with  variable  coeffients  has  been  developed. 

Presented  in  Figs.  4  and  5  are  the  three  dimensional  plots  for  the  results  of 
e^(y,  n)  and  e,/y,  m),  where  h  =  2  is  assumed.  Using  the  relations  [23, 
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24]  derived  from  Faraday’s  law 
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the  corresponding  components  of  magnetic  field  perturbation  By,  B^  and  |B  |  are 

de, 

also  obtained  and  presented  in  figures  6-8,  where  =  1  x  10  is  assumed  at 
the  magnetopause.  Using  the  same  procedure  for  the  stability  analysis  for  m  =  20 
case,  the  optimum  Ai/  =  0.028  is  found.  It  is  then  incorporated  in  the  numerical 
scheme  in  solving  (1)  and  (2)  for  m  =  20.  The  results  of  and  together 
with  B^  and  B^  obtained  through  (10) -(12)  are  presented  in  Figs.  9-13. 
The  results  show  that  the  field  line  resonance  (toroidal)  mode  (£y,B^,BM) 
can  be  excited  inside  the  magnetosphere  through  the  coupling  with  the  global  cavity 
(poloida!)  mode  (e^,  By)  initiated  by  the  perturbation  appeared  on  the  outer 
boundary  of  the  magnetosphere  (i.e.,  magnetopause).  The  field  line  resonance  mode 
appears  to  be  a  localized  oscillation  as  expected.  For  m  =  20,  the  peak  of  the 
oscillation  is  located  along  v  =  0.23.  It  moves  outward  to  v  =  0.13  when  m  is 
increased  to  100.  It  is  noted  that,  for  m  =  100  case,  the  frequency  considered  is  too 
low  for  the  globe  mode  to  propagate  in  the  magnetosphere.  Nevertheless,  through  the 
tunneling  effect  the  field  line  resonance  mode  can  still  be  excited  by  coupling  to  the 


excited  by  coupling  to  the  decay  part  of  the  globe  mode  oscillation. 

VII.  Summary  and  Conclusion 

A  numerical  algorithm  dealing  with  the  mixed-type  coupled  partial  differential 
equations  with  variable  coefficients  has  been  developed.  The  stability  of  the 
numerical  scheme  has  been  examined  and  solutions  for  the  hydromagnetic  wave 

equations  are  obtained  on  non-staggered  grids  by  this  algorithm.  The  main 

contribution  of  this  algorithm  is  to  offer  an  efficient  way  to  examine  the  stability 
conditions  first  by  the  Fourier  method  and  then  check  them  by  the  matrix  method. 

Since  these  two  methods  are  different,  one  can  be  sure  of  the  stability  conditions. 

Moreover,  this  algorithm  leads  a  way  for  a  parametric  study  of  the  stability 
conditions.  It  also  becomes  an  advantage  in  the  sense  that  the  optimum  parameters 
which  can  increase  the  accuracy  of  the  numerical  solution  can  be  determined  together 
with  the  stability  analyses.  Using  the  developed  numerical  scheme,  the  problem  of 
coupling  and  propagation  of  hydromagnetic  waves  in  the  realistic  dipole  model  of  the 
magnetosphere  can  be  studied.  An  example  has  been  considered.  The  results  are 
presented  in  Figs.  4-13. 
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Figure  Captions 

Fig.  1.  Maximum  eigenvalue  evaluated  by  the  Fourier  method  along  the  v  axis  at 
A#i  =  0.004,  Ai/  =  0.01334  for  different  6. 

Fig.  2.  Maximum  eigenvalue  evaluated  by  the  Fourier  method  along  the  v  axis  at 
A/i  =  0.004,  9  =  1  for  different  Av. 

Fig.  3.  Maximum  eigenvalue  evaluated  by  the  Matrix  method  along  the  v  axis  at 
9  =  1,  Am  =  0.004,  Ai/  =  0.01334. 

Fig.  4.  Normalized  electric  field  of  poloidal  mode  (m  =  100). 

Fig.  5.  Normalized  electric  field  of  toroidal  mode  (m  =  100). 

Fig.  6.  Normalized  magnetic  field  of  poloidal  mode  (radial  component)  (m  = 
100). 

Fig.  7.  Normalized  magnetic  field  of  toloidal  mode  (N-S  component)  (m  =  100). 

Fig.  8.  Normalized  magnetic  field  of  toroidal  mode  (E-W  component)  (m  = 

100). 

Fig.  9.  Normalized  electric  field  of  poloidal  mode  (m  =  20). 

Fig.  10.  Normalized  electric  field  of  toroidal  mode  (m  =20). 

Fig.  II.  Normalized  magnetic  field  of  poroidal  mode  (radial  component)  (m  =  20). 

Fig.  12.  Normalized  magnetic  field  of  toroidal  mode  (N-S  component)  (m  =  20). 

Fig.  13.  Normalized  magnetic  field  of  toroidal  mode  (E-W  component)  (m  = 

20). 
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genvalue  evaluated  by  the  Fourier  method  along  the  v  axis  at 
6  =  1  for  different  Ay. 


-0.00286 


genvalue  evaluated  by  the  Matrix  method  along  the  v  axis  at 
0.004.  =  0.01334. 


Normalized  magnetic  field  of  toroidal  mode  (N-S  component) 
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Abstract.  The  temporal  evolution  of  Langmuir 
waves  excited  by  hlgh-power,  high-frequency  (HF) 
radio  waves  In  the  Ionosphere  Is  studied  theoreti¬ 
cally.  This  study  Is  motivated  by  past  observations 
made  with  the  450  MHz  radar  at  Areclbo  Observa¬ 
tory  In  Puerto  Rico.  Two  kinds  of  nonlinear  damp¬ 
ing  to  the  parametric  decay  instability  are  considered 
In  the  derivation  of  the  rate  equation  for  the  spectral 
Intensity  of  enhanced  Langmuir  waves.  They  are 
Langmuir  wave  cascading  caused  by  nonlinear  Lan¬ 
dau  damping  and  cross-field  electron  diffusion.  The 
first  damping  process  leads  to  the  saturation  of  Indi¬ 
vidual  unstable  Langmuir  wave.  The  second  pro¬ 
cess,  which  results  from  the  incoherent  scattering  of 
electron  orbits  by  the  total  excited  Langmuir  waves, 
yields  anomalous  damping  that  applies  to  each  Lang¬ 
muir  wave.  Consequently,  Langmuir  waves  with 
smaller  growth  rates  will  be  suppressed  by  those  with 
larger  growth  rates.  Such  a  mode  competition  pro¬ 
cess  may  cause  the  overshoot  of  the  HF-enhanced 
plasma  line  observed  with  the  Areclbo  430  MHz 
radar.  Favorable  agreement  is  obtained  between 
theory  and  the  Areclbo  observations. 

Introduction 

The  HF-enhanced  plasma  lines  observed  in  the 
Areclbo  heating  experiments  have  been  generally 
attributed  to  the  parametric  decay  Instability, 
whereby  the  plasma  waves  and  heavily  damped  lon- 
acoustlc  wave  are  parametrically  excited  by  the 
Injected  HF  heater  waves  [e.g.  Fejer  and  Kuo,  1973; 
Perkins,  et  al.,  1974].  The  excited  plasma  waves 
then  enhance  the  backscatter  spectrum  of  the  Are¬ 
clbo  430  MHz  radar  at  frequencies  near  430  MHz 
±  fHF  where  fHF,  Is  the  frequency  of  the  HF  heater 
wave.  The  radar  returns  at  these  two  sidebands  are 
often  referred  to  as  HF-enhanced  plasma  lines 
(HFPL,).  One  of  the  most  reproducible  phenomena 
associated  with  the  HF-enhanced  plasma  lines  In  the 
F  region  over  Areclbo  is  the  so  called  main  plasma 
line  overshoot  [Showen  and  Kim,  1978;  Showen  and 
Behnk-,  1978].  A  clear  demonstration  of  the 
overshoot  phenomenon  Is  presented  In  the  Figure  2 
of  Showen  and  Behnke.  However,  the  growth  and 
decay  periods  of  the  overshoot  presented  In  that  Fig¬ 
ure  are  rather  long  In  comparison  with  the  recent 
observations  at  Areclbo.  It  Is  shown  In  the  recent 
study  of  temporal  evolution  of  HFPLS  that  the  HFPL 
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signal  exhibited  an  Initial  growth  lasting  for  a  few 
tens  of  milliseconds  before  reaching  a  maximum 
value,  and  then  greatly  reduced  its  strength  In  the 
order  of  a  second.  The  typical  time  scales  for  the 
overshoot  process  of  HFPL,  at  Areclbo  are  summar¬ 
ized  In  Table  I.  These  values  are  based  upon  the 
recent  data  presented  by  Duncan  and  Sheerln  [1985] 
and  DJuth  et  al.  [1986],  and  the  other  unpublished 
data.  Based  upon  the  experimental  observation  of 
Coster  et  al.,  [1985]  at  Areclbo,  It  seems  unlikely 
that  the  anomalous  absorption  by  HF-!nduced  shortr 
scale  (1-lOm)  Irregularities  can  account  for  this 
overshoot  phenomenon.  In  addition,  recent  efforts 
made  at  Areclbo  to  carefully  measure  the  differential 
absorption  of  the  lonospherlcally  reflected  HF  wave 
have  found  an  upper  bound  for  the  anomalous 
absorption  of  the  HF  waves  by  the  Induced  shortr 
scale  Irregularities  that  Is  no  more  than  10% power  at 
the  heater  frequency  or  5.1  MHz  [J.A.  Fejer,  F.T. 
DJuth,  and  M.P.  Sulzer,  private  communications, 
1988].  Given  the  fact  that  strong  plasma  line 
overshoot  was  present  during  the  absorption  meas¬ 
urements,  this  data  supports  the  Idea  that  anomalous 
absorption  does  not  play  a  key  role  In  the  plasma¬ 
line  overshoot  phenomenon. 

In  the  present  work  we  propose  a  mechanism 
that  Is  based  on  the  mode  suppression  process  Intro¬ 
duced  by  the  resonance  broadening  effect  to  account 
for  the  plasmarline  overshoot  phenomenon.  It  is 
expected  that  the  presence  of  the  parametric  decay 
Instabilities  gives  rise  to  perturbations  on  the  phase 
space  orbits  of  electrons  (charged  particles  In  gen¬ 
eral).  The  result  of  the  Incoherent  scattering  of  elec¬ 
tron  orbits  by  the  total  excited  plasma  waves  can 
lead  to  electron  diffusion  In  the  velocity  space  along 
the  magnetic  fleid  together  with  the  cross-field 
diffusion  In  the  coordinate  space.  On  one  hand, 
these  diffusion  processes  appear  as  an  enhanced 
viscosity  to  the  electron  motion.  But,  on  the  other 
hand,  they  broaden  the  resonance  Interaction 
between  the  electrons  and  waves.  It  has  been  sug¬ 
gested  that  the  resonance  broadening  effect  can  con¬ 
tribute  to  the  saturation  of  the  parametric  Instability 
[Weinstock  and  Bezzerldes,  1973]  mainly  through 
the  velocity  diffusion  process,  and  further  lead  to  the 
mode  suppression  [Dupree,  1988;  Dum  and  Dupree, 
1970]  effectively  through  the  cross-field  spatial 
diffusion  process.  The  mode  suppression  process  can 
be  understood  from  the  fact  that  the  anomalous 
damping  Introduced  by  resonance  broadening  effect 
to  each  unstable  plasma  line  Is  proportional  to  the 
total  spectral  Intensity  of  the  excited  plasma  waves. 
In  other  words,  the  enhanced  electron  diffusion  Is 
determined  by  the  total  spectral  Intensity.  The 
increase  of  the  anomalous  damping  on  the  slowly 
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TABLE  I.  The  Arecibo  430  MHz  HF-Enhanced 
Plasma  Line  (HFPL)  In  the  F  Region 


Range 
of  Values 

Typical 

Value 

Total  HFPL 

Growth  Period 

2-50  ms 

20  ms 

Overshoot 

Time  Scale 

0.1  - 

10.0  s 

1  s 

HFPL  Overshoot 

Ratio 

10 

-  100 

50 

growing  plasma  line  can  result  from  the  fast  growth 
of  other  lines  In  the  same  region.  It  will  be  shown 
later  that  the  linear  growth  rate  of  each  plasma  line 
depends  on  the  wave  propagation  angle  8,  (the 
geomagnetic  field)  and  the  wavenumber  k  of  the 
plasma  line.  Consequently,  the  spectral  lines  having 
smaller  growth  rates  will  be  suppressed  when  their 
linear  growth  rates  are  exceeded  by  the  anomalous 
damping  introduced  by  the  modes  with  larger  growth 
rates.  A  physical  model  based  on  the  aforesaid  pro¬ 
cess  will  be  presented  In  the  subsequent  sections  to 
interpret  the  plasma-line  overshoot  phenomenon. 

Theoretical  Model 

We  consider  the  parametric  decay  Instability  near 
the  reflection  height  of  the  HF  heater  wave  in  the 
ionospheric  F-region.  If  an  HF  heater  wave  with  the 
ordlnaj-y  polarization  Is  transmitted,  the  wave  electric 
field  Ep  =  2Ep  can  be  assumed  to  be  parallel  to  the 
geomagnetic  field  B0  =  2B0  In  the  region  near  Its 
reflection  height.  Electron  plasma  waves  and  Ion 
acoustic  waves  are  then  excited  parametrically  by  the 
imposed  pump  wave  as  the  pump  field  Intensity 
exceeds  the  threshold  of  the  Instability.  The  linear 
dispersion  relation  of  this  Instability  for  a  dipole 
pump  wave  Ep  =  E0e~ IWot  +  c.c.  Is  given  by  [Nlshi- 
kawa,  1968;  Kuo  and  Cheo,  1978], 

4tusu;(i(T +^/2)  (T+P,/2)  =  k2(e2E2/mem,)  (l) 

where  uis  and  uig  represent  the  frequencies  of  the  Ion 
acoustic  mode  and  electron  plasma  mode,  respec¬ 
tively;  the  frequency  matching  condition 
cjs  +  u.'q  =  u>0  Is  satisfied;  7  Is  the  linear  growth  rate, 
Ug  Is  the  electron-ion  collision  frequency; 
I/,  =  (7r/2)l/2a/s  ( Te/T,)3/2e  T,/2Tr 3,2  represents  the 
ion  Landau  damping  rate,  and  kt  =  kcos#  for  the 
excited  spectral  line  with  k  and  a  propagation  angle  9 
with  respect  to  the  earth  magnetic  field. 

The  linear  growth  rate  7^  of  the  spectral  line  k 
Is  then  round  to  be 

TkL(t=  -j-(  +  [;»E0/2(m,m,u;oug,/2] 

cos#  =  «k1,,2cos#-  0O-  0xk  (2) 

where  a  =  eE0/2(mem,u;0cs)l/2,  30  =  izg/4. 

Pi  =  -J-(tr/2)1/2cs  (Te/T,)3/2e  T'/2Tr3/2,  ws  =  kcs  and 

utg~u)0  have  been  used,  and  the  maximum  linear 
growth  rate  is  given  by  7^0max  =  (o'2/-lf?1)  cos2#  /?,, 


for  k  =  k0  =  a2cos2#/4/?2.  It  is  seen  that  the  linear 
growth  rate  of  the  spectral  line  is  a  function  of  the 
propagation  angle  9  and  the  wave  number  k  and  the 
growth  rate  is  maximized  for  each  k  at  9  =  0.  How¬ 
ever,  the  instability  can  also  be  excited  at  a  large 
propagation  angle  If  the  pump  wave  Is  Intense 
enough.  For  a  pump  field  exceeding  the  threshold 
field,  the  maximum  angle  for  the  unstable  |k  |  modes 
is  given  by 

#ok  =  cos"  1[(^0  +  ^1k)/ok1/2]  (3) 

where  0O+Pl  k<  «k1/2  is  required  for  the  instabil¬ 
ity. 

We  now  consider  two  possible  nonlinear  damping 
mechanisms  of  the  Instability.  One  Is  due  to  the  cas¬ 
cading  or  the  plasma  lines  and  the  other  one  Is  the 
cross-field  diffusion  damping  caused  by  the  reso¬ 
nance  broadening  efTect.  The  cascading  process 
occurs  when  the  amplitude  of  the  excited  plasma  line 
exceeds  the  threshold  of  parametric  decay  process. 
The  daughter  plasma  line  then  takes  energy  away 
from  the  mother  plasma  line  and  thus  slows  down 
the  growth  of  the  mother  plasma  line.  If  this  decay 
process  does  not  drain  much  energy  from  the 
mother  plasma  line,  cascading  process  will  proceed 
further  to  excite  more  decay  lines  until  the  total 
energy  drained  by  the  decay  lines  Is  large  enough  to 
stop  the  growth  of  the  mother  line.  Based  on  the 
energy  balance  relation,  the  anomalous  damping  rate 
t/k*0  Introduced  by  the  cascading  process  on  the  origi¬ 
nal  decay  Instability  is  approximately  equal  to  the 
collisionless  linear  growth  rate  of  the  daughter  decay 
mode  having  the  mother  decay  mode  as  the  pump. 
It  is 

pk*#=  v/k  eEk<,/2(mem,w0cs)1/2  (4) 

where  Ekff  is  the  field  amplitude  of  the  mother  decay 
plasma  line  propagating  at  the  angle  8  and  having  the 
wavenumber  k. 

In  addition  to  the  cascading  process  that  can  only 
lead  to  the  saturation  or  the  instability,  a  second 
non-linear  effect  which  originates  in  the  broadening 
of  wave-electron  resonances  by  the  random  motion 
of  electrons  in  the  excited  turbulent  field  Is  also  con¬ 
sidered  in  the  study  of  the  HFPL  overshoot 
phenomenon.  Secular  changes  of  the  guiding  center 
positions,  cyclotron  radii,  and  phase  angles  of  elec¬ 
trons  by  the  turbulent  field  give  rise  to  resonance 
broadening  and  cross-field  diffusion  similar  to  those 
produced  by  collisions!  scattering.  It  has  been  shown 
by  Dum  and  Dupree  ]1970]  that  the  field-dependent 
broadening  can  be  expressed  in  terms  of  diffusion 
coefficients  which  broaden  the  resonance  function  In 
a  similar  way  as  the  collisions!  damping  does.  Con¬ 
sequently,  an  anomalous  damping  to  the  instability  Is 
also  introduced  by  the  resonance  broadening  process. 
In  the  magnetized  plasma,  the  cross-field  diffusion 
will  be  the  major  response  of  electrons  to  the 
Incoherent  scattering  by  waves.  The  diffusion  damp¬ 
ing  rate  to  the  k  line  is  thus  obtained  as 

d„(k)  =  kj/Dj/kL  (5) 

where  u  =  £(c2(ZxEk,)  (2xEk.)/B2)  T~^  7  • 

f.  (<  +  T  k2,o) 

7k0  =  k/iPocos#',  and  O  2«^2  ls  turned.  Replac¬ 
ing  the  summation  by  integral,  we  can  rewrite  (5)  as 
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Fig.  1.  Temporal  evolution  of  the  spectral  inten¬ 
sity  of  the  *10°  plasma  lines;  the  parameters  2E0 
=  0..'>v/m,  p,,  =  400/sec,  f0  =  5MHz  and  Cs  =  3.5xl03 
m/sec  are  used. 

d„(k)  =  [ac2k2sin20/B2tu2I 
/  /  2?rk'2  vie7  sin3  0'costf'  e£„,  dk'  d9'.  (6) 


~@o  =  ~  {bsin30ocos',0o/[a  sin20o  -I-  bISosins0ocoss0o 
+  ( 2-  5sin20o/cos20o)(  a#ocos0o  -  3fo)}  } 

ej,  t) 

J  [agcosff  -  fig  -  2dtf-  t}elg!2\  I9sin2  9  cos  6  d  9 

O 

(12) 

where  b  =  a6c2/32/?fB2w 2  and  a  =  a2/0l. 

So  far  we  have  derived  two  coupled  equations 
(9)  and  (12)  for  the  temporal  evolution  of  the  spec¬ 
tral  intensity  of  HFPLS  oriented  at  the  angle  6 .  They 
have  been  solved  numerically.  For  describing  the 
overshoot  phenomenon  observed  In  the  Arecibo 
heating  experiments,  we  let  9  =  40  *  and  the  numeri¬ 
cal  results  of  (9)  and  (12)  are  presented  in  Figures  1 
and  2,  which  show  the  functional  dependence  of  I40- 
and  90(t)  on  time  t,  respectively.  As  shown  In  Fig¬ 
ure  1,  the  spectral  Intensity  I40-  grows  Initially  at  the 
linear  growth  rate.  It  overshoots  for  very  high  level 
to  form  a  peak  and  then  decreases  at  a  much  slower 
rate  subsequently.  It  eventually  saturates  at  a  much 
lower  level  than  Its  peak  value.  An  overshoot 
phenomenon  has  been  clearly  demonstrated  by  the 
proposed  theory. 


Fronj  (2),  (4)  and  (6),  the  nonlinear  growth  rate  of 
the  k  line  can  be  defined  by 


7k9  =  7ke~  <*e(V-  (7) 

The  rate  equation  for  describing  the  nonlinear  evolu¬ 
tion  of  Eka  Is  then  given  by 

~TT^k 9  ~  (8) 


Analysts 

To  simplify  the  analysis,  we  only  analyze  the 
evolution  of  the  total  spectral  intensity  1$  of  plasma 
lines  propagating  at  the  same  angle  8  with  respect  to 
the  geomagnetic  field,  where  I#  =  4trsin0/  dk  k2Ek$ 
is  so  defined.  In  terms  of  I9,  Equation  (8)  becomes 
approximately 


Qg  COS0  -  0e~  2dfl  -  rjell'2 


Io 


(0) 


where  ne  =  2av/k7>  3$  =  2{30  +3,  k0).  Ve 


=  e/k01/2  (4n’sin^mem,u;0cs)1 


and 


d„  = 


ac2 k® /2sln20/2B2u>2J  f  ls,sin20'cost)’dO'; 

0 

k0  =  a2cos20/4;?2,  (10) 


we  perform  the  average  in  k  by  assuming  the  Ek29  is 
peaked  Initially  at  k0  with  maximum  linear  growth 
rate;  dQ  is  the  maximum  averaged  spectral  angle  and 
defined  by 


'Wos0„  -  39o  -  2d#0  =  0  (ll) 

With  the  aid  of  (9)  and  (10),  taking  the  time 
derivative  of  (11)  leads  to 


Summary 

A  theoretical  model  for  the  temporal  evolution 
of  the  HF-enhanced  plasma  line  (HFPL)  Is 
developed.  We  suggest  that  mode  suppression 
caused  by  the  resonance  broadening  effect  may  be 
the  physical  process  responsible  for  the  overshoot  of 
HFPL  detected  by  the  Arecibo  430  MHz  radar.  In 
other  words,  the  presence  of  other  decay  lines  hav¬ 
ing  larger  growth  rates,  which  cannot  be  detected  by 
the  Arecibo  radar,  gives  rise  to  enhanced  viscosity  to 
electron  motion  and,  thus.  Increases  the  threshold 
for  exciting  the  detected  plasma  line.  We  also 
include  the  nonlinear  effect  Introduced  by  the  cas¬ 
cading  process  in  the  formulation  theory.  This  effect 
results  in  the  saturation  of  the  parametric  decay  Ins¬ 
tability,  and  hence,  Is  responsible  for  the  steady  state 
level  of  the  HFPL  after  the  overshoot  process. 

The  difference  between  the  linear  growth  rate  of 
the  parametric  decay  instability  and  the  anomalous 
damping  Introduced  the  aforementioned  two 
processes  defines  the  nonlinear  growth  rate  of  the 


Fig.  2.  Evolution  of  the  maximum  unstable  angle 
in  time. 
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HFPL.  It  Is  then  used  for  establishing  the  nonlinear 
growth  rate  equation  (8)  of  the  HFPL,  from  which 
the  rate  equation  (8)  for  the  total  spectral  Intensity  le 
of  plasma  lines  propagating  at  the  same  angle  8  Is 
derived.  Although  Equation  (0)  does  not  provide 
the  exact  amount  of  the  temporal  development  of 
HFPL,  It  greatly  simplifies  the  analysis  of  the  prob¬ 
lem  and  maintains  the  Important  features  of  the 
phenomenon  as  Illustrated  In  Figure  1.  A  second 
rate  equation  (12)  for  the  averaged  maximum  spec¬ 
tral  angle  60  is  also  obtained  for  studying  the  tem¬ 
poral  evolution  of  the  HFPL  self-conslstently.  The 
numerical  results  of  (8)  and  (12)  for  9  =  40 '  have 
clearly  demonstrated  the  overshoot  phenomenon 
similar  to  that  observed  by  Showen  and  Behnke 
(1878J  in  the  Arecibo  Ionospheric  heating  experi¬ 
ments.  However,  regarding  the  overshoot  time 
scales,  our  calculated  values  are  more  consistent  with 
the  typical  results  of  recent  observations  at  Arecibo 
as  shown  by  Table  I. 
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1-our  invariants  of  the  ray  trajectory  are  found  for  a  ray  propagating  in  a  horizontally  stratified 
ionosphere  under  the  density  perturbation  of  HF  wave-induced  field-aligned  irregularities  The  reflec¬ 
tion  height  of  the  ray  can  then  be  determined  with  the  aid  of  those  invariants.  The  results  show  that  the 
reflection  height  of  the  ray  varies  drastically  (namely,  strong  spread  F  echoes)  in  the  presence  of 
irregularities  that  polarize  in  the  magnetic  meridian  plane  By  contrast,  the  reflection  height  is  not 
affected  (namely,  no  spread  F  echoes)  by  those  irregularities  that  polarize  in  the  direction  perpendicular 
to  the  meridian  plane.  Spread  F  is  quite  insensitive  to  the  magnetic  dip  angle  80  in  the  region  from  20 
to  70  .  The  dependence  of  spread  F  on  the  scale  length  of  the  irregularity  has  also  been  examined  for 
the  case  0„  =  50\  It  is  found  that  spread  F  is  not  caused  by  irregularity  with  scale  length  less  than 
about  100  m. 


1  INTRODUCTION 

Spread  F  that  refers  to  diffuse  echoes  on  an  iono- 
gram  from  the  ionospheric  F  region  was  observed  in 
the  ionospheric  heating  experiments  conducted  at 
Boulder,  Colorado  [ Utlaut ,  1971].  This  ionospheric 
phenomenon  is  generally  believed  to  be  caused  by 
the  excitation  of  large-scale  (a  few  hundreds  of  meters 
to  kilometers)  field-aligned  ionospheric  irregularities. 
However,  spread  F  is  a  rare  phenomenon  in  the  ex¬ 
periments  at  Arecibo,  Puerto  Rico  [Showed  and  Kim. 
1978]  and  it  has  not  been  observed  at  Tromso, 
Norway  since  the  new  European  heating  facility  has 
been  in  operation  [Stuhhe  et  a!.,  1982].  Evidences, 
such  as  the  radio  star  scintillations  [Frey  et  ai,  1984] 
and  the  sanning  radar  incoherent  backscatter  process 
[Duncan  and  Behnke.  1978]  indicate  that  large-scale 
ionospheric  irregularities  have  been  excited  by  HF 
heater  waves  at  Tromso  and  Arecibo.  Hence  a  lack 
of  spread  F  echoes  does  not  imply  the  absence  of 
heater  wave-induced  ionospheric  irregularities.  It 
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may  be  due  to  the  difference  in  the  polarization  di¬ 
rections  of  the  HF  wave-induced  irregularities  where 
the  polarization  direction  of  the  irregularity  is  de¬ 
fined  to  be  the  direction  of  its  wave  vector.  In  gener¬ 
al,  field-aligned  irregularities  can  have  two  indepen¬ 
dent  polarization  directions.  One  lies  in  the  meridian 
plane  and  the  other  one  is  in  the  direction  perpen¬ 
dicular  to  the  meridian  plane.  The  theoretical  results 
of  filamentation  instability  in  magnetoplasmas  [Kuo 
and  Schmidt.  1983]  also  show  that  the  irregularities 
excited  by  the  o  mode  pump  and  by  the  x  mode 
pump  have  different  polarization  directions.  The  ir¬ 
regularities  excited  by  the  o  mode  pump  are  field- 
aligned  and  are  polarized  in  the  direction  perpen¬ 
dicular  to  the  meridian  plane.  By  contrast,  the  irregu¬ 
larities  excited  by  the  x  mode  pump  are  polarized  in 
the  meridian  plane  and  are,  in  general,  not  field- 
aligned.  However,  the  field-aligned  nature  of  the  ir¬ 
regularities  may  be  established  to  reduce  the  diffu¬ 
sion  damping  along  the  magnetic  field. 

Since  spread  F  is  generally  used  as  an  indicator  of 
large-scale  irregularities,  and  yet  it  is  known  that 
sometimes  such  irregularities  do  not  produce  spread 
F.  Hence  a  study  leading  to  better  understanding  of 
the  spread  F  phenomenon  and  its  relationship  with 
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Fig-  1.  Coordinates  used  in  the  analysis. 


density  irregularities  is  thought  to  be  important  in 
this  aspect.  In  the  present  work  a  theoretical  model 
for  artificial  spread  F  echoes  is  developed,  from 
which  the  relationship  between  the  spread  F  echoes 
and  C  HF  wave-induced  irregularities  is  studied. 
The  primary  purpose  of  this  study  is  to  determine  the 
effects  of  the  irregularity  polarizations,  scale  length, 
and  the  magnetic  dipangle  on  the  spread  F  echo.  In 
section  2,  the  ray  trajectory  equations  are  analyzed 
to  study  the  wave  propagation  in  a  horizontally  stra¬ 
tified  ionosphere  under  the  density  perturbation  of 
the  HF  wave-induced  field-aligned  irregularities. 
Four  invariants  of  the  ray  trajectory  are  found.  Main 
results  are  illustrated  in  the  figures.  Conclusions  are 
finally  drawn  in  section  3. 

2  MODEL  AND  ANALYSIS 

We  consider  a  horizontally  stratified  ionosphere 
having  a  scale  length  JL.  Thus  the  unperturbed  elec¬ 
tron  density  is  represented  by  «0(_k)  =  !  +  x/L), 

where  x  is  the  vertical  coordinate  and  n0  =  n0( 0)  is 
the  electron  density  at  the  reference  plane  v  =  0  lo¬ 
cated  at  height  H.  When  field-aligned  irregularities 
are  present,  the  background  electron  density  distri¬ 
bution  is  perturbed.  Thus  the  total  density  is  the  sum 
of  the  unperturbed  density  and  the  fluctuating  den¬ 
sities  of  the  irregularities.  In  general,  field-aligned  ir¬ 
regularities  can  have  two  independent  polarization 
directions.  One  lies  in  the  magnetic  meridian  plane 
and  the  other  one  is  in  the  direction  perpendicular  to 
th?  meridian  plane.  For  those  irregularities  which  arc 
field-aligned  in  the  magnetic  meridian  plane,  the  total 
density  perturbation  is  simply  expressed  in  a  form  of 
sin  A; ,(.v  cos  fl0  +  z  sin  ft0).  where  2n  k,  =  is 
the  average  scale  length  of  the  irregularities.  0o  is  the 
magnetic  dip  angle,  and  z  is  the  horizontal  coordi¬ 
nate  as  shown  in  Figure  1.  The  density  perturbation 
associated  with  the  irregularities  which  are  field 
aligned  in  the  direction  perpendicular  to  the  mag¬ 


netic  meridian  plane  is  modeled  as  dn2  sin  {k2y  +  <£), 
where  /.2  =  2zr  k2  is  the  averaged  scale  length,  y  is  the 
coordinate  perpendicular  to  the  meridian  plane,  and 
<t>  is  the  arbitrary  phase  angle.  In  this  model,  the  total 
electron  density  distribution,  n(.v)  is  then  composed 
of  three  components,  namely. 

tt(x.  z)  —  n0(l  +  v  Li  +  dn,  sin  A,(.v  cos  tl0  +  :  sin  0a} 

+  Sn2  sin  (k2  y  +  <t>)  (I) 

For  simplicity,  the  effect  of  the  geomagnetic  field 
on  wave  propagation  will  be  neglected  in  the  follow¬ 
ing  analyses.  The  dispersion  relation  for  the  wave 
propagation  is  thus  given  bv 

ID 2  =  ID +  k2c2  (2) 

where  id  and  k  are  the  wave  frequency  and  the  wave 
number,  respectively,  and 

ai2pt  =  <o£,0[l  +  '  F  +  (<$n(/n0)  sin  fc,(x  cos  0o  4-  z  sin  0o) 

+  (dn2  n0)  sin  (k2  y  +  &)]  (3| 

Since,  in  reality,  the  ionosonde  antennas  radiate 
radio  wave  beams  with  finite  cross  sections  and 
spread  angles,  so  it  is  more  appropriate  to  consider 
the  transmitted  pulses  to  be  composed  of  rays  whose 
initial  wave  vectors  and  the  horizontal  locations  on 
the  reference  plane  x  =  0  are  all  different.  The  trajec¬ 
tory  of  each  ray  is  then  governed  by  the  following  set 
of  characteristic  equations: 


d 

—  x  =  k.  C2!ID 

dl 

(4) 

d  , 

—  l  =  k..  i'Jid 
dt  ’ 

(5) 

d 

—  z  —  k.c2 do 
dl 

(6| 

-(<n,,0  2d>L)[l  4-  (dn,  n„Uk,L  cos  (t„| 


cos  k  ,(x  cos  ()0  +  z  sin  (ln)] 


(7| 


—  kt  =  ~U’>2e„  2oH^n2  n„)k2  cos  (L,  y  +  4> I  (Si 


-tan  0„{iD2r„  2o>l.<  ,Sn,  A , cos  ft,,) 


cos  fc,(x  cos  (tn  +  z  sin  ()0)  (9) 

where  (51  and  (8)  form  one  set  of  coupled  equations, 
and  (4).  (6),  (7)  and  (9)  form  the  other  set  of  coupled 
equations,  but  these  two  sets  of  equations  do  not 
couple  to  each  other.  However,  those  coupled  equa- 
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tions  can  be  separated  by  the  invariants  of  the  trajec¬ 
tory  derived  as  follows.  We  first  take  the  ratio  of  (5) 
and  (8),  which  yields 


dy_ 
dk . 


2k,  t 2 


(10) 


",  “>Uo(&ni/no)ki  cos  (k2y  +  <t>) 

Integrating  this  equation  gives  the  first  invariant  of 
the  trajectory 

(Sn2/n0)  sin  (k2y  +  <p)  +  k2c2/<ol, o  =  const  in  time 

=  (Srt2/n0)  sin  (k2  y0  +  <j>)  +  k20  c2/a>l,0  (11) 

where  the  subscript  0  represents  the  initial  value. 

We  next  substitute  (9)  into  (7)  to  yield 


—  (k*  ~  cot  e„k,)  =  -to1  Ill'll 
<ti 


(12) 


Integrating  (12)  from  the  initial  height  x  =  0  of  the 
ray  to  the  new  height  x  at  time  t,  leads  to  the  second 
invariant  relation 

k,  —  cot  0O  k,  +  (a>l.0/2a>L)t  =  const  in  time 

~  kgo  cot  dpljo  (13) 

We  now  construct  an  equation  by  first  multiplying 
(4)  and  (6)  by  cos  0O  and  sin  0O  respectively,  and  then 
summing  them  up.  The  result  is 

d 

—  ( x  cos  0o  +  z  sin  60)  =  (c2/eoXk,  cos  0o  +  k,  sin  0O) 


(14) 

Two  additional  invariants  will  be  derived  as  follows. 
From  the  ratio  of  (14)  and  (7),  we  obtain 

d(x  cos  90  +  1  sin  0o)  2 c2L 

dk,  <0 

k,  cos  0o  +  k,  sin  0a 

I  +  (Snt/n0Xk,L  cos  0O)  cos  M-x  cos  0O  +  z  sin  0o) 

(15) 

With  the  aid  of  (13),  (14)  and  (4),  (15)  is  integrated  to 
obtain  the  third  invariant 

.x  +  (tan  280/2)z  +  (<5n,/it0)L(cos2  0O  cos  20o) 

sin  k  ,(x  cos  60  +  z  sin  60) 

+  (c2L/ro£,0)<fcj  +  tan  20okxk,) 

=  const  in  time  =  (tan  20o/2)zo  +  (6 njn0) 

U cos2  0o,/cos  20o)  sin  (fc , r0  sin  0O) 

+  [c1Liio2p.0Xk2x0  +  tan  2 80kx0kr0)  (16) 

The  fourth  invariant  is  obtained  from  integrating  the 


resultant  equation  defined  by  the  ratio  of  (14)  and  (9). 
It  is 

(tan  20o/2)z  +  (<5n,/n0)L(sin2  0o/cos  20o) 

■  sin  k,(x  cos  60  +  z  sin  0O) 

+  (c2L/w^M-kf  +  tan  20ofc,fcI) 

=  const  in  time  =  (tan  280/2)z0  +  (Snl/nQ) 

■  Lfsin2  B0/ cos  20o)  sin  (k,z0  sin  0O) 

+  (c2L/(jj£,0X  -  k20  +  tan  280  kx0  k,0)  ( 1 7) 

However,  (16)  and  (17)  can  be  rearranged  to  reduce 
to  much  simpler  forms.  We  first  take  the  difference  of 

(16)  and  (17)  to  yield 

x  +  (3n,/n0)L  sin  fc,(x  cos  80  +  z  sin  0o )  +  (k2c2 L/co^o) 

=  (8n,/n0)L  sin  (k,z0  sin  0O)  +  kgC2L/(o^o  (18) 

where  k0  =  k(0).  We  next  multiply  (16)  by  sin2  0O  and 

(17)  by  cos2  90  and  then  take  their  difference.  The 
result  is 

x  -  z  cot  0O  +  (cIL/a>^r0Ykx  -  k,  cot  0O)2 

=  -z0  cot  0„  +  (c2L/t<jJ,oKkro  ~  k,o  cot  0O)2  (19) 


Fig.  2.  Example  of  the  appropriate  trajectories  of  the  rays  de¬ 
tectable  to  the  ionosonde.  The  parameters  used  tn  the  calculation 
are  9„  =  50',  L  =  50  km,  a,  =  1  km,  kof2/<uJ»o  =  *•  and  1  ” 
An,  h0  =  0.1. 
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Fig.  3.  The  functional  dependence  of  the  reflection  height  spread  of  a  diagnostic  beam  on  the  magnetic  dip 

angle  0o  for  nn  =  (a)  0.005,  ( b )  0.01,  (r)  0.05,  and  (<f)  0. 1 . 


We  have  shown  that  the  trajectory  of  each  ray  is 
governed  by  four  invariant  relations  (11),  (13),  (18) 
and  (19)  that  are  derived  from  two  sets  of  coupled 
differential  equations  (4H9).  Hence  the  temporal  ev¬ 
olution  of  any  one  of  the  four  variables  x,  kx  and 
k..  and  either  one  of  the  two  variables  y  and  k,. .  can 
be  determined,  with  the  aid  of  these  invariants,  by 
the  corresponding  rate  equations,  namely,  equations 
(4)  (9).  The  elapsed  time  for  each  ray  traveling  from 
the  reference  height  to  the  reflection  height  can  in 
principle  be  determined  by  integrating  (7)  from  kx  = 
kt0  to  kx  =  0,  with  the  prescribed  initial  conditions: 
x(0)  =  0,  -(0)  =  :0  and  k.(0)  =  k.0.  Since  the  density 
perturbation  tin2  produced  by  the  irregularities  ori¬ 
ented  in  the  direction  perpendicular  to  the  magnetic 


meridian  plane  does  not  appear  in  (7),  the  reflection 
height  of  each  ray  is  thus  not  expected  to  be  modified 
by  the  presence  of  such  oriented  irregularities.  Physi¬ 
cally.  it  is  realized  that  a  ray  is  reflected  (i.e.,  from 
upgoing  to  downgoing)  when  kx  changed  the  sign. 
We  thus  expect  that  the  reflection  height  of  the  ray  is 
likely  determined  by  those  parameters  which  can 
affect  the  value  of  kx .  Apparently,  equation  (7)  shows 
that  (>n2  is  not  the  one  among  those  parameters.  This 
result  that  kx  is  independent  of  Sn2  is  also  consistent 
with  the  Snells'  law  in  the  limit  of  the  absence  of  the 
v  variation  in  the  plasma  density.  However,  the  den¬ 
sity  perturbation  i in,  produced  by  the  irregularities 
that  are  oriented  within  the  meridian  plane  can  affect 
the  reflection  height  of  the  ray  as  seen  in  (7).  One 
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Fig.  4.  Axm„  versus  x  for  0o  =  (a)  50  .  (h)  68  ,  and  (<■)  78  .  where 
x  =  dn ,  n0 . 


way  to  show  this  effect  is  to  integrate  (7)  numerically. 
This  equation  can  be  integrated  with  the  aid  of  the 
invariant  relations  (13),  (18)  and  (19).  If  we  assume 
that  each  ray  of  the  antenna  radiation  beam  has  a 
auaight  trajectory  from  the  ionosonde  to  the  refer¬ 
ence  plane  as  shown  in  Figure  1.  then  the  initial 
conditions  for  ray  trajectory  starting  from  the  refer¬ 
ence  plane  can  be  defined  as  x  =  x0  ,  kx0  —  k0  H/(H2 
+  Zq)1  \  and  kx0  =  k0z0/(H 2  +  Zq)1  2,  where  z0  is  a 
free  parameter  characterizing  the  ray.  By  varying  z0 
within  the  cross  sectional  range  of  th«,  beam,  a  set  of 
trajectories  for  the  rays  of  the  beam  can  be  deter¬ 
mined.  However,  not  all  the  rays  can  be  reflected 


back  to  the  ionosonde.  Only  the  rays  which  can 
return  to  the  ionosonde  are  responsible  for  the  ob¬ 
served  diffuse  echoes  on  the  ionogram,  and  hence  the 
intensity  of  spread  F  is  proportional  to  the  maximum 
difference  of  the  reflection  heights  of  those  rays.  The 
boundary  conditions  at  x  =  0  for  a  downgoing  ray 
which  can  propagate  straightly  from  the  reference 
plane  back  to  the  ionosonde  are  defined  to  be  kx  = 
-  k0H/(H2  +  :2f)il2  and  kx  =  -k0zf/(H2  +  z})112, 
where  zf  is  the  z  coordinate  of  the  ray  on  the  refer¬ 
ence  plane.  Those  conditions  can  be  satisfied  only  by 
rays  with  appropriate  z0.  In  other  words,  the  con¬ 
ditions  at  x  =  0  for  the  upgoing  rays  which  can  be 
received  later  by  the  ionosonde  can  not  be  prede¬ 
fined.  Instead,  z0  and  zf  have  to  be  determined  simul¬ 
taneously. 

Presented  in  Figure  2  are  the  two  distinctive  tra¬ 
jectories  for  two  different  groups  of  rays  in  the  beam, 
where  the  background  parameters  used  in  the  calcu¬ 
lation  are  0O  =  50°,  L  =  50  km,  x,  =  1  km,  k%c2/ 
to^o  =  1  and  &njn0  =  a  =  0.1.  The  rays  which 
follow  these  trajectories  can  return  to  the  ionosonde. 
These  two  trajectories  are  found  to  depend  on  the 
value  of  a.  They  approach  each  other  and  then  merge 
into  a  single  trajectory  coincided  with  the  x  axis  as 
the  value  of  a  decreases  to  zero.  It  is  therefore  shown 
that  if  the  ionosonde  transmitting  beam  is  modelled 
by  many  rays  having  different  initial  locations  on  a 
reference  plane,  the  spread  F  echoes  can  indeed  be 
interpreted  to  be  caused  by  the  difference  in  the  re¬ 
flection  heights  (or  the  time  delays)  of  the  returned 
signals.  The  virtual  height  spread  is  thus  proportion¬ 
al  to  the  maximum  difference  Axm.„  of  these  reflec¬ 
tion  heights.  This  spread  disappears  when  the  irregu¬ 
larities  oriented  in  the  meridian  plane  are  absent.  It 
should  be  noted  that  an  ionosonde,  like  any  receiver, 
has  a  limited  sensitivity  and  requires  a  certain  mini¬ 
mum  “number  of  rays”  from  any  height  interval  to  be 
recognized  as  a  signal.  Therefore,  the  reflection 
height  of  this  signal  should  be  defined  to  be  the 
average  value  over  those  of  rays.  In  view  of  this, 
Axm„  defined  in  this  study  may  not  be  a  fully  sensi¬ 
ble  measure  of  the  spread  F  strength.  Nevertheless,  it 
indeed  provides  a  qualitative  indication. 

We  now  employ  the  similar  parameters  to  analyze 
the  system  of  (7),  (13),  (18),  and  (19):  L  =  50  km, 
:  i  —  1  km  and  k$c2 /a>^0  =  1 ;  z0  is  allowed  to  vary 
for  determining  the  trajectories  of  the  rays.  Only 
taken  into  account  are  the  trajectories  of  the  rays 
detectable  to  the  ionosonde  for  determining  the  max¬ 
imum  difference  A.v.m„  from  the  reflection  heights  ot 
these  rays.  Axm,„  versus  the  magnetic  dip  angle  0o  as 
calculated  for  a(  =  Snt/n0 )  =  0.005,  0.01,  0.05  and  0.1 
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Fig.  5.  The  dependence  of  spread  F  on  the  scale  length  x,  of  the 
irregularity,  where  60  =  50  and  x  =  0.05. 

and  presented  in  Figures  3a~3d.  Shown  in  Figures 
4a~4c  are  Axm„  versus  the  irregularity  intensity  a  for 
0O  =  50°,  68°,  and  78°  that  correspond  to  magnetic 
dip  angles  at  Arecibo,  Boulder  and  Tromso,  respec¬ 
tively.  The  dependence  of  spread  F  on  the  scale 
length  of  the  irregularity  has  also  been  examined  for 
the  case  80  =  50°  and  a  =  0.05.  The  result  is  present¬ 
ed  in  Figure  5  showing  that  spread  F  can  only  be 
caused  by  irregularities  with  scale  lengths  larger  than 
about  106  m. 

3  DISCUSSION  AND  CONCLUSIONS 

A  theoretical  model  of  spread  F  echoes  has  been 
developed  foi  determining  the  effects  of  the  magnetic 
dip  angle,  irregularity  intensity,  and  the  polarization 
direction  and  the  scale  length  of  irregularity  on  the 
occurrence  frequency  and  the  intensity  of  spread  F. 
We  first  show  that  irregularities  polarized  perpen¬ 
dicular  to  the  meridian  plane  do  not  cause  spread  F. 
Therefore,  only  irregularities  polarized  within  the 
meridian  plane  have  been  considered  in  the  analyses 
for  evaluating  the  aforementioned  effects. 

As  shown  in  Figures  3a-3d,  spread  F  becomes 
prominent  as  the  magnetic  dip  angle  exceeds  about 
5°.  The  reflection  height  spread  Axm„  first  increases 
rapidly  to  a  maximum  at  about  60  =  8°  (except  for 
very  low  irregularity  intensity,  e.g.,  a  =  0.001,  where 
the  peak  appears  at  0O  s  20°)  and  then  decreases 
monotonically  to  zero  at  90  =  90°.  However,  Axm„  is 
quite  flat  in  the  region  from  20°  to  70°  (30°  to  80°  for 
a  =  0.001).  Furthermore,  Axm„  increases  monotoni¬ 
cally  with  the  irregularity  intensity  as  seen  in  Figures 
4a,  4 b  and  4c  for  90  =  50°,  68°,  and  78°  correspond¬ 
ing,  respectively,  to  the  magnetic  dip  angles  at  Are¬ 
cibo,  Boulder,  and  Tromso.  While  a  near  linear  de¬ 
pendence  is  found  in  the  0O  =  50°  and  68°  cases, 
Axma,  starts  to  approach  a  maximum  at  a  =  0.02  for 
60  =  78c.  The  dependence  of  spread  F  on  the  scale 


length  A,  of  irregularities  has  been  examined  for 
0O  =  50°  and  a  =  0.05.  The  result  shown  in  Figure  5 
indicates  that  spread  F  appears  only  when  A,  >  106 
m.  The  reflection  height  spread  remains  a  fairly  con¬ 
stant  value  from  106  m  to  about  500  m.  This  value 
then  increases  monotonically  by  35%  at  A,  =  1  km 
which  is  the  largest  scale  length  employed  in  the 
analyses.  Our  model  shows  that  spread  F  is  quite 
insensitive  to  the  magnetic  dip  angle.  As  long  as  the 
irregularity  with  the  polarization  within  the  meridian 
plane  exists,  significant  spread  F  should  be  observed 
over  a  wide  range  of  latitude  including  Arecibo, 
Boulder,  and  even  Tromsa.  In  other  words,  the  ap¬ 
pearance  of  spread  F  should  not  depend  upon  the 
locations  of  the  ionosondes.  We,  therefore,  conclude 
that  spread  F  echoes  on  the  ionograms  are  caused  by 
the  irregularities  with  polarization  directions  within 
the  meridian  plane  and  the  scale  lengths  greater  than 
about  100  m.  By  contrast,  irregularities  with  polar¬ 
ization  directions  perpendicular  to  the  meridian 
plane  do  not  give  rise  to  spread  F  echoes. 

The  different  occurrence  frequencies  of  artificial 
spread  F  noticed  at  Arecibo,  Boulder,  and  Tromso 
are  most  probably  due  to  the  excitation  of  different 
types  of  irregularities.  However,  the  large  geomagnet¬ 
ic  dip  angle  at  Tromso  also  plays  an  important  role 
(see  Figure  3),  which  works  to  substantially  reduce 
the  spread  F  strength  at  Tromso.  The  o  mode  and  x 
mode  pump  waves  transmitted  from  the  Boulder 
heating  facilities  cannot  be  separated  as  easily  as 
those  from  the  Arecibo  or  the  Tromso  facilities.  This 
fact  can  be  evidenced  by  the  measurements  of  HF 
wave-induced  short-scale  irregularities  at  Boulder 
[ Fialer ,  1974]  showing  that  x  mode  can  still  excite 
short-scale  irregularities  though  not  as  efficiently  as  o 
mode  wave.  The  heater  wave  transmitted  into  the 
ionosphere  from  Tromso.  Norway,  is  a  left  (x  mode) 
or  right  (o  mode)  hand  circularly  polarized  wave 
propagating  along  the  geomagnetic  field.  Therefore 
the  generated  irregularities  (by  either  o  mode  or  x 
mode)  at  Tromso  are  expected  to  be  magnetic  field- 
aligned  [Kuo  and  Schmidt ,  1983],  Further,  if  the 
plasma  inhomogeneity  in  the  vertical  direction  is 
taken  into  account,  the  filamentation  instability 
would  be  excited  preferentially  in  the  direction  per¬ 
pendicular  to  the  meridian  plane  (i.e.,  the  y  direction 
of  the  coordinate  system  used  in  this  work)  in  order 
to  have  symmetric  high  frequency  sidebands.  This, 
together  with  the  effect  of  geomagnetic  dip  angle, 
may  explain  why  spread  F  has  never  been  oh«*rved 
since  the  Tromso  facilities  were  operated  a  tew  years 
ago  We  also  note  that  when  the  heater  is  operated  in 
o  mode  at  Arecibo.  no  spread  F  or  change  of  reflec- 
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tion  heights  can  be  seen  (L.  M.  Duncan,  private  com¬ 
munication,  1984).  It  then  implies,  based  on  the  pro¬ 
posed  model  of  spread  F  mechanism,  that  the  o  mode 
heater  excited  irregularities  must  be  oriented  in  the  y 
direction.  This  again  agrees  with  the  theoretical  pre¬ 
diction  of  the  generation  of  large-scale  irregularities 
by  the  filamentation  instability  [ Kuo  and  Schmidt, 
1983]  and  supports  the  proposed  model  of  spread  F 
mechanism  self-consistently. 
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We  elaborate  upon  the  physical  mechanism  of  thermal  filamentation  instability  of  radio  waves  whose 
frequencies  can  be  as  low  as  in  the  VLF  band  and  as  high  as  in  the  SHF  band.  This  instability  can 
excite  large-scale  magnetic  and  plasma  density  fluctuations  simultaneously  in  the  ionosphere  and 
magnetosphere.  We  comment  on  relevant  experiments  in  terms  of  this  instability  and  other  mecha¬ 
nisms. 


1  INTRODUCTION 

Manifest  ionospheric  disturbances  can  be  pro¬ 
duced  by  powerful  HF  radio  waves  in  ionospheric 
heating  experiments,  such  as  fluctuations  in  iono¬ 
spheric  density,  plasma  temperature,  and  the  earth’s 
magnetic  field  (see,  for  example,  Gurevich  [1978], 
Fejer  [1979],  Stubbe  et  al.  [1982],  and  Lee  and  Kuo 
[1983a]).  Unexpectedly  large  perturbations  in  the 
earth’s  magnetic  field  { — 10.8  nT)  were  observed  in 
Tromso  experiments  with  the  E1SCAT  HF  heating 
facilities  [ Stubbe  and  Kopka,  1981],  and  we  predicted 
the  coexistence  of  large-scale  ionospheric  irregu¬ 
larities  [Kuo  and  Lee,  1983].  Whistler  waves  at  fre¬ 
quencies  close  to,  but  less  than,  the  local  electron 
gyrofrequency  are  expected  to  cause  both  plasma 
density  irregularities  and  geomagnetic  field  fluctu¬ 
ations  in  the  ionosphere  if  MF  signals  are  transmit¬ 
ted,  or  in  the  magnetosphere  if  VLF  signals  are  in¬ 
jected  [ Lee  and  Kuo ,  1984a].  Microwave  transmis¬ 
sions  at  2.45  GHz  with  central  maximum  power  den¬ 
sity  of  230  W/m2  from  the  conceptualized  Solar 
Power  Satellite  are  also  expected  to  perturb  the 
earth's  magnetic  field  significantly  in  the  ionosphere 
along  the  microwave  beam  [ Lee  and  Kuo,  1984b]. 
Our  theoretical  analyses  have  shown  that  geomag¬ 
netic  field  fluctuations  can  be  excited  simultaneously 
with  large-scale  field-aligned  ionospheric  density  by 
powerful  radio  waves.  The  frequencies  of  radio  waves 
Copynght  1985  by  the  American  Geophysical  1  Inion. 
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can  be  as  low  as  in  the  VLF  band  and  as  high  as  in 
the  SHF  band. 

Thermal  filamentation  instability  excitea  by  radio 
waves  in  the  ionosphere  and  magnetosphere  is  re¬ 
sponsible  for  the  simultaneous  perturbation  in  iono¬ 
spheric  density  SN  and  geomagnetic  field  6B.  A  non- 
oscillatory  (i.e.,  purely  growing)  mode  is  associated 
with  3N  and  bB,  whose  magnitudes  are  interrelated 
uniquely  by  an  equation  characteristic  of  the  thermal 
filamentation  instability.  Since  the  filamentation 
wave  vector  is  perpendicular  to  the  incident  wave 
vector,  the  radio  wave  tends  to  be  “self-focused”  with 
this  instability.  One  purpose  of  this  paper  is  to  elab¬ 
orate  upon  the  physical  mechanism  of  thermal  fila¬ 
mentation  instability  as  presented  in  section  2,  that 
leads  to  features  distinctively  different  from  those  of 
thermal  self-focusing  instabilities  [e.g.,  Perkins  and 
Valeo,  1974;  Cragin  et  al.,  1977;  Perkins  and  Gold¬ 
man ,  1981].  In  addition,  we  wish  to  comment  on 
some  relevant  experiments  in  section  3  in  terms  of 
thermal  filamentation  instability  and  other  mecha¬ 
nisms.  Conclusions  are  finally  drawn  in  section  4. 

2  PHYSICAL  MECHANISM  OF  THERMAL 
FILAMENTATION  INSTABILITY 

Radio  waves,  if  intense  enough,  are  able  to  excite 
thermal  filamentation  instability  in  the  ionosphere 
and  magnetosphere,  that  produces  sideband  and 
purely  growing  modes  with  wave  vectors  perpendicu¬ 
lar  to  the  propagation  direction  of  radio  waves.  Mag¬ 
netic  field-aligned  nonoscillatory  modes  (i.e.,  fila- 
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Fig.  L  The  physical  process  leading  to  the  simultaneous  exci¬ 
tation  of  large-scale  earth's  magnetic  field  perturbations  <5B  and 
plasma  density  irregularities  SN  by  the  thermal  filamentation  in¬ 
stability  of  powerful  radio  waves. 

ments  of  magnetostatic  and  ionospheric  density  fluc¬ 
tuations)  are  formed  along  the  radio  wave  paths. 
They  are  caused  by  a  nonlinear  thermal  effect  due  to 
radio  wave-plasma  interactions.  We  elucidate  its 
physical  mechanism  as  follows  and  diagram  it  in 
Figure  1. 

The  wave  fields,  that  include  the  incident  and  side¬ 
band  fields,  interact  with  plasmas  and  consequently 
yields  two  types  of  nonlinear  force:  the  pondcr- 
omotive  force  (or  generally  termed  nonlinear  Lorentz 
force)  and  the  thermal  pressure  force.  In  generating 
field-aligned  nonoscillatory  modes,  the  thermal  pres¬ 
sure  force  turns  out  to  be  larger  than  the  ponder- 
omotive  force  by  a  factor  of  (A/2nr,)2  for  modes  with 
scale  lengths  A  less  than  nr,  (2 M/m)112  (~  15  m  in  the 
ionosphere)  and  by  a  factor  of  (M/m)  otherwise  [Kuo 
et  al.,  1983;  Lee  and  Kuo ,  19836],  where  r,  and  (M/m) 
are  the  electron  gyroradius  and  the  ratio  of  ion  to 
electron  masses,  respectively.  Ions  cannot  be  ef¬ 
fectively  heated  because  of  large  cross-field  heat  diffu¬ 
sion  loss.  Thermal  pressure  force  can  be  built  up  only 
in  the  electron  gas.  The  dotted  arrow  in  Figure  1 
denotes  such  a  force  fT  =  -x{d/dx)(N0ST,)  experi¬ 
enced  by  electrons  across  the  earth’s  magnetic  field 
(B0  =  zB0). 

Electron  density  fluctuations  (SN)  that  are  induced 
consequently  by  the  thermal  pressure  force  give  rise 
to  a  self-consistent  field  (<5E  =  xSE )  associated  with 
the  excitation  of  the  purely  growing  mode.  While  the 
net  force  acting  on  electrons  is  (fr  —  eN0  <5E),  that  on 
singly  charged  ions  is  eN0SE.  These  wave-induced 
forces  combined  with  the  orthogonal  earth’s  mag¬ 
netic  field  result  in  cross-field  drift  motions  of  elec¬ 
trons  and  ions  at  the  velocities  of  (<5E  —  fr/ejV0) 
x  Bo  'Bo  an<3  (<5E  x  B0/Bj),  respectively.  Therefore. 
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Fig  2.  The  threshold  condition  and  a  positive  feedback  loop  for 
the  thermal  filamentation  instability. 


a  relative  motion  between  electrons  and  ions  at  the 
velocity  of  (fr/eN0)  x  B 0/Bl  results  and  produces  a 
net  current  flowing  along  the  direction  of  the  y  axis. 
This  wave-induced  current  (and  vector  potential,  <5A) 
varies  sinusoidally  along  the  x  axis.  It  can  be  visual¬ 
ized  as  a  combination  of  many  sheetlike  currents  on 
the  y-z  plane  with  a  thickness  of  half  the  scale  length 
A,  flowing  in  opposite  directions  alternately.  Or, 
equivalently,  pairs  of  dipole  current  form  within  the 
radio  wave-heated  ionospheric  region.  Earth’s  mag- 
entic  field  fluctuations  <5B  that  result  from  the  radio 
wave-induced  current  are  field  aligned.  They  point  to 
the  same  direction  as  the  geomagnetic  field’s,  B0, 
because  of  V  x  SA  =  AB  =  zSB. 

The  physical  mechanism  of  producing  both  6B  and 
SN  is  summarized  in  Figure  2  with  the  aid  of  block 
diagrams,  showing  how  a  positive  feedback  loop  can 
be  formed  for  the  thermal  filamentation  instability. 
In  brief,  electron  temperature  perturbation  ST, 
caused  by  radio  wave  heating  produces  a  cross-field 
thermal  pressure  force  fr  as  the  key  nonlinearity  of 
the  thermal  filamentation  instability.  This  nonlinear 
force  induces  earth’s  magnetic  field  perturbation  6B 
and  ionospheric  density  irregularities  SN  that,  subse¬ 
quently,  cause  nonuniform  electron  heating  and, 
then,  perturbs  the  electron  temperature  in  turn.  It 
can  be  expected  from  the  above  delineated  physical 
picture  that  SN  x  SB  because  both  of  them  are  in¬ 
duced  in  proportion  by  fr .  Such  a  relation  has  been, 
indeed,  found  in  our  detailed  formulation  of  the 
theory  with  the  following  form  [equation  (2  )  of  Lee 
and  Kuo,  19846] : 


SN 
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where  the  parameters,  v^v,),  0,(0,),  wpe ,  c,  >•,  a  •_(/.„), 
and  ,V0(£0)  have  their  conventional  meanings  of  ef¬ 
fective  electron  (ion)  collision  frequency,  electron 
(ion)  gyrofrequency,  electron  plasma  frequency,  the 
speed  of  light  in  vacuum,  the  growth  rate  of  the  in¬ 
stability,  the  perpendicular  (parallel)  scale  length  of 
the  instability  and  the  background  ionospheric  den¬ 
sity  (earth’s  magnetic  field  strength),  respectively.  The 
density  perturbation,  SN/N0 ,  is  proportional  linearly 
to  the  magnetic  perturbation,  SB/B0 ,  with  a  ratio 
greater  than  unity.  If  this  ratio  is  close  to  unity, 
SB/B0  comparable  to  SN/N0  can  be  excited  by  the 
thermal  filamentation  instability.  Otherwise,  the 
earth’s  magnetic  field  perturbation  is  negligibly  small 
as  compared  with  the  ionospheric  density  pertur¬ 
bation. 

One  might,  however,  be  skeptical  about  the  possi¬ 
bility  of  producing  SB/B0  ^  SN/N0  based  on  the  fol¬ 
lowing  argument.  If  SB/B0  ~  SN/N0 ,  then  both 
6B/B0  and  SN/N0  must  be  extremely  small  compared 
with  STr/T0  according  to  the  “MHD  equilibrium” 
condition: 


M  SB  P  ST, 
N o  '  B0  '  2  To 


(2) 


where  0  is  the  ratio  of  kinetic  to  magnetic  pressure, 
namely,  N0  T0/(Bq  .'8  it)  ~  10  5  for  typical  iono¬ 
spheric  parameters.  It  is  apparent  from  (2)  that,  in 
"MHD  equilibrium,”  STe;T0  (~2  x  I02)  is  absurdly 
unrealistic  even  for  0.1%  of  <j/V  /V0  or  SB/B0.  The 
argument  that  SB/B0  must  be  extremely  small  under 
the  “MHD  equilibrium”  condition  is  certainly  true. 
This  fact  can  be  also  seen  from  (1),  which  requires 
SB  J30  =  0  when  y  =  0,  corresponding  to  the  "MHD 
equilibrium"  condition  either  before  the  onset  or 
after  the  saturation  of  the  thermal  filamentation  in¬ 
stability.  During  the  linear  stage  of  the  instability 
(i.e..  the  growth  rate  (y)  =  a  nonzero  constant),  SB/B0 
and  6S;  ,V0  are  related  by  ( I )  and  significant  geomag¬ 
netic  field  fluctuations  can  be  produced  by  powerful 
radio  waves.  In  other  words,  the  geomagnetic  field 
perturbation  caused  by  powerful  radio  waves  via  the 
thermal  filamentation  instability  is  a  transient  phe¬ 
nomenon,  whose  duration  may  be  estimated  from  the 
growth  rate  of  the  instability.  It  is  reasonable  to 
define  this  duration  to  be  the  period  for  achieving  the 
seven  e-folds  of  magnitude  above  the  thermal  fluctu¬ 
ation  level,  namely,  7y  '. 

The  parallel  wave  number  k  is  zero  for  ideal  field- 
aligned  modes.  Nevertheless,  a  small  (compared  to 
k  }  but  finite  k  can  be  introduced  by  inhomogeneity 


effects  associated  with  either  the  background  iono¬ 
spheric  density  or  the  finite  cross  section  of  the  radio 
wave  beam.  The  inhomogeneity  effects  have  been 
reasonably  ignored  in  our  previous  studies  of  VLF, 
MF,  and  HF  cases.  This  is  because  the  instability 
becomes  effective  in  the  ionosphere  and  mag¬ 
netosphere  for  modes  with  scale  lengths  less  than  the 
scale  sizes  of  the  background  plasma  density  gradi¬ 
ents  and  the  linear  dimensions  of  wave  beam  cross 
sections  by,  at  least,  2  orders  of  magnitude.  By  con¬ 
trast,  the  inhomogeneity  effect  imposed  by  the 
narrow  microwave  beam  cannot  be  neglected  in  the 
case  of  the  Solar  Power  Satellite.  The  parallel  scale 
length  a,i  ,  which  is  of  the  order  of  the  microwave 
beam  size  (~  10  km),  renders  the  factor  enclosed  by 
the  brackets  in  (1)  to  be  much  greater  than  unity  in 
the  ionospheric  F  region.  Physically,  large  heat  con¬ 
duction  loss  along  the  geomagnetic  field  inhibits  the 
thermal  filamentation  instability  from  producing  sig¬ 
nificant  earth's  magnetic  field  perturbation  in  the  F 
region.  Parallel  heat  conduction  loss  is,  however, 
drastically  reduced  by  the  large  electron-neutral  col¬ 
lisions  in  the  E  region.  Therefore,  the  thermal  fila¬ 
mentation  instability  can  operate  effectively  in  the  £ 
region  but  not  the  F  region  in  the  solar  power  satel¬ 
lite  case. 

Parallel  heat  conduction  loss  may  significantly  en¬ 
hance  the  threshold  of  the  instability  whose  full  ex¬ 
pression  is  given  by  [equation  (6  )  of  Lee  and  Kuo. 
1984b] 


where  e.,h,  a±(a  ).  m(M),  and  V,  are  the  threshold  wave 
field  intensity,  the  perpendicular  (parallel)  scale 
length,  the  electron  (ion)  mass,  and  the  electron  ther¬ 
mal  velocity,  respectively;  the  factor,/,  can  be  a  func¬ 
tion  of  a>0,  iope,  Qe.  and  /.L  depending  upon  the 
frequency  (co0)  of  the  radio  waves.  The  thermal  fila¬ 
mentation  instability  requires  a  threshold  basically 
because  the  sideband  and  purely  growing  modes  are 
not  normal  modes  but  nonlinearly  driven  modes  in 
space  plasmas.  The  threshold  level  (c,h)  is  determined 
in  (3).  representing  the  energy  balance  between  the 
driving  thermal  source  of  the  instability  (i.e.,  the  col- 
lisional  dissipation  of  radio  wave  energy)  and  the 
damping  processes  denoted  by  the  three  terms: 
2 (m'MI,  (2n  Vt  a1  Qp)2,  and  (2nV,//.v  v,)2.  These  three 
damping  terms  correspond  to  the  collisional  damp¬ 
ing  of  the  thermal  source,  the  cross-field  heat  conduc¬ 
tion  loss,  and  the  parallel  heat  conduction  loss,  re- 
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spectively.  The  first  term  dominates  for  the  instability 
with  >  (2M/m)lllnVJtlt  (~!5  m  in  the  iono¬ 
sphere)  and  (An/AJ  >  (Ot/vt).  But,  parallel  heat  con¬ 
duction  loss  becomes  the  dominant  damping  process 
if  A||/Amfp  >  ?r(2M/m)1/2,  where  AMFP  is  the  electron 
mean  free  path. 

In  general,  the  instability  with  large-scale  lengths 
requires  lower  thresholds  since  1  £,* |  oc  1 .  The  fa¬ 
vorable  excitation  of  geomagnetic  field  fluctuations 
by  large-scale  modes  can  be  seen  in  (1),  which  for 
highly  field  aligned  modes,  can  be  reduced  to 


6N 


y  \A<v/  J\b0 


(4) 


where  the  growth  rate  y  increases  generally  with  scale 
lengths.  Smaller  ratios  of  SN/N0  to  SB/B0  are  there¬ 
fore  associated  with  larger-scale  modes,  namely,  sig¬ 
nificant  geomagnetic  field  fluctuations  can  be  excited 
mainly  through  large-scale  thermal  filamentation  in¬ 
stability. 

It  is  interesting  to  note  that  the  threshold  con¬ 
dition  given  by  (3)  is  independent  of  the  electron  col¬ 
lision  frequency  v,  for  the  excitation  of  “large-scale, 
highly  field  aligned”  modes,  that  is. 


ffii  f4nV<  ( 
me  t.L  \M ) 


(5) 


where  the  factor  /  is  not  a  function  of  v, .  This  is 
understandable  because  both  the  generation  rate  (the 
left-hand  side  of  (5))  and  the  damping  rate  of  the 
thermal  source  are  proportional  to  v, .  More  specifi¬ 
cally.  electron  collisions  with  ions  or  neutrals  convert 
the  radio  wave  energy  into  the  thermal  source  of  the 
filamentation  instability;  meanwhile,  the  induced 
electron  temperature  perturbations  are  also  deter¬ 
mined  by  the  collisional  loss  of  wave-generated  heat. 
It  should  be  stressed  that  electron  collisions  do  not 
weaken  the  purely  growing  modes  (i.e.,  the  mag¬ 
netostatic  and  ionospheric  density  fluctuation)  be¬ 
cause  electrons  and  ions  move  together  in  the  wave 
vector  direction  of  the  purely  growing  mode. 

The  expressions  of  the  growth  rate  have  been 
found  to  be  generally  proportional  to  v, .  This  offers 
another  evidence  that  electron  collisions  enhance  the 
instability.  The  principal  reason  is  that  large  cross¬ 
field  mobility  of  charged  particles  is  conducive  to  the 
establishment  of  collective  oscillations,  i.e.,  the  field- 
aligned  purely  growing  fluid  modes.  The  cross-field 
mobility  is  greatly  enhanced  by  electron  collisions. 
Another  interesting  thing  to  note  is  that  the  relation¬ 
ship  between  SN  and  SB  shown  in  (4)  turns  out  to  be 
not  a  function  of  v,  because  y  x  vr . 


The  generation  of  large-scale  ionospheric  density 
irregularities  by  powerful  radio  waves  have  been  in¬ 
vestigated,  for  example,  by  Perkins  and  Valeo  [1974] 
and  Perkins  and  Goldman  [1981]  via  the  thermal  self- 
focusing  instability,  and  by  Fejer  [1973],  Cragin  and 
Fejer  [1974],  and  Cragin  et  al.  [1977]  via  the  ther¬ 
mal  stimulated  Brillouin  backscattering  instability. 
These  two  instabilities  have  their  own  distinctive  fea¬ 
tures  in,  for  instance,  the  threshold  conditions  dis¬ 
cussed  in  next  section,  though  they  are  believed  to  be 
physically  equivalent.  The  simultaneous  excitation  of 
large-scale  earth's  magnetic  field  fluctuations  and 
ionospheric  density  irregularities  is  the  important 
discovery  in  our  theory  of  thermal  filamentation  in¬ 
stability.  In  addition,  the  geomagnetic  field  effect  on 
the  high-frequency  sideband  mode  that  was  ignored 
in  the  previous  theories  has  been  taken  into  account 
in  our  formulation  of  the  instability. 

The  inclusion  of  this  effect  exhibits  its  important 
contribution  to  the  driving  sources  of  sideband 
modes  (especially  as  <o0  <  fie)  that  satisfy  the  follow¬ 
ing  Fourier  transform  wave  equation 

+  k2  -  -yjl  -  k^zz  -  k2 ix  -  i k„  —  lx:  -  ixtjt, 

An  «„ 

—  I  — - —  (J5  +■  J \  +  J:«l  (6) 

c 

where  I  is  the  unit  dyadic;  f  and  x  denote  the  unit 
vectors  taken  along  the  direction  of  the  geomagnetic 
field  (B0  =  “80)  and  the  filamentation  wave  vector 
(k  =  xk).  For  simplicity,  the  radio  wave  has  been  as¬ 
sumed  to  propagate  along  the  geomagnetic  field.  The 
w  uve-induced  current  density  includes  three  parts: 
Js=—eN0Xl,  Js=eNX0,  and  Jw  =  e.V0VM, 
where 


V,= 


1  r>n 
m  (im2  -  ft2) 


j  =  0  or  I  (7) 

are  the  electron  quiver  velocities  responding  to  the 
incident  wave  field  (e0)  and  the  sideband  field  (e,). 
respectively; 
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is  the  electron  velocity  perturbation  caused  by  the 
V0  x  <SB  Lorentz  force.  Js  is  the  linear  response  of 
electrons  to  the  sideband  field  e,.  whereas  J„  and  JM 
are  the  nonlinear  beating  currents  stemming  from  the 
nonlinear  coupling  between  the  purely  growing 
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modes  (AN  or  <>B)  and  the  incident  wave  field  (£0).  On 
one  hand,  the  sideband  mode  is  produced  by  the 
parametric  interaction  of  the  incident  radio  wave  and 
the  purely  growing  modes.  On  the  other  hand,  the 
purely  growing  modes  are  driven  by  the  thermal 
pressure  force  ( —  V{N0ATe))  that  is  caused  by  the 
ohmic  dissipation  of  the  incident  and  sideband 
modes  at  the  rate  of  Qe  =  2if,V0m(VJ  •  V,  +  V0  • 

vr>. 

J.  COMMENTS  ON  RELEVANT  EXPERIMENTS 

The  thermal  filamentation  instability  of  VLF 
waves  (whistlers)  can  be  excited  in  the  magnetosphere 
rather  than  the  ionosphere.  This  is  because  the  insta¬ 
bility  can  only  be  excited  in  the  frequency  nnge 
n,  2  <  o>„  <  Q,  by  whistlers  and  Q,/27t  is  about  1.4 
MHz  in  the  ionosphere  and  13.65  KHz  in  the  mag¬ 
netosphere  at  L  =  4.  Although  VLF  wave  injection 
experiments  have  been  performed  actively,  for  exam¬ 
ple.  at  Siplc,  Antarctica,  the  ground-based  transmit¬ 
ters  are  usually  operated  in  a  pulsed-wave  mode  with 
duration  of  a  few  seconds.  These  VLF  wave  pulses 
have  been  used  to  study,  primarily,  coherent  wave- 
particle  interactions  in  the  magnetosphere  such  as 
the  wave  amplification,  the  triggering  of  wave  emis¬ 
sions,  the  induced  particle  precipitation,  etc.  (see,  for 
example,  Hellinell  [1983]).  However,  the  excitation 
of  thermal  filamentation  instability  in  the  mag¬ 
netosphere  requires  the  continuous  operation  of  VLF 
transmitters  for  a  few  minutes  because  of  the  small 
electron  collision  frequency  ( ~  0. 1  Hz)  and  then  the 
small  growth  rates  (10  3  to  10  2  Hz).  Experiments 
can  be  planned  with  the  Siple  transmitter  or  more 
powerful  U  S.  Navy  communication  transmitters  op¬ 
erated  in  C'W  modes  and  also  with  spacecrafts,  for 
instance,  the  ISEE  I  to  monitor  the  VLF  wave- 
induced  magnetospheric  disturbances. 

The  thermal  filamentation  instability  can  be  excit¬ 
ed  in  the  ionosphere  by  radio  waves  whose  fre¬ 
quencies  exceed  half  the  electron  gyrofrequcncy 
(-0.7  MHz).  The  ionospheric  heating  facilities  lo¬ 
cated  at  Arecibo  and  Tromso  currently  operate  in 
the  HF  band  at  lowest  frequencies  of  the  order  of  3 
MHz  (J.  A.  Fejer,  private  communication,  1984). 
Those  in  U  S  S  R,  are  also  primarily  HF  transmitters 
except  the  one  near  Moscow  that  operates  at  a  fre¬ 
quency  (-1.35  MHz)  near  the  electron  gy- 
rofrequency  [ Gurevich  and  Migulin.  1982].  Unex¬ 
pectedly  large  geomagnetic  field  perturbation  of  the 
order  of  10  nT  had  been  caused  by  the  Tromso  trans¬ 
mitter  [Siubbe  and  Kopka,  1981;  Stuhbe  et  al..  1982] 
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and  we  explain  it  in  terms  of  the  thermal  fila¬ 
mentation  instability  [Kuo  and  Lee,  1983].  This  in¬ 
stability  needs  several  minutes  to  be  developed  under 
the  experimental  conditions  at  Tromso.  This  charac¬ 
teristic  growth  time  is  consistent  with  the  observa¬ 
tions  that  geomagnetic  perturbation  increased  regu¬ 
larly  with  the  operation  time  of  radio  waves  in  the 
range  from  10  to  360  s  [ Stubbe  et  al.,  1982J.  Another 
interpretation  of  magnetic  fluctuations  is  that  they 
are  caused  by  perturbations  in  the  steady  currents 
(i.e.,  the  auroral  electrojet)  flowing  in  the  ionosphere 
caused  by  perturbations  in  the  ionospheric  conduc¬ 
tivity  [ Fejer  and  Krenzien,  1982,  and  references 
therein].  One  would  expect  from  this  mechanism  that 
increased  magnetic  field  fluctuations  are  associated 
with  more  intense  electrojet  current  for  future  experi¬ 
mental  verification.  The  source(s)  of  magnetic  fluctu¬ 
ations  can  be  unambiguously  identified  by  in  situ 
measurements  performed  by  rockets  or  spacecrafts  in 
the  modified  ionosphere.  While  the  auroral  electrojet 
as  a  line  source  is  located  at  an  altitude  of  about  100 
km  above  the  earth's  crest,  the  radio  wave-induced 
dc  current  appears  in  a  higher  and  broader  heated 
ionospheric  region  below  the  reflection  height  of  the 
incident  heater  wave.  A  unique  feature  of  thermal 
filamentation  instability  is  that  the  simultaneously 
excited  magnetic  and  density  fluctuations  have  iden¬ 
tical  characteristic  scale  lengths  and  growth  rates. 
Moreover,  the  strengths  of  these  two  types  of  fluctu¬ 
ations  are  interrelated  specifically  by  (1).  The  in  situ 
measurements  would  provide  an  excellent  diagnosis 
of  the  heated  ionosphere  to  test  the  theories  in  these 
aspects. 

The  theories  of  Perkins  and  Valeo  [1974]  and 
Cragin  et  al.  [1977]  have  been  examined  by  Farley  et 
al.  [1983]  and  Frey  and  Duncan  [1984]  as  candidate 
mechanisms  of  producing  large  scale  irregularities  in 
overdense  ionospheric  heating  experiments  at  Are¬ 
cibo.  The  thresholds  of  Perkins  and  Valeo  theory  are 
significantly  higher  than  those  of  Cragin  et  al.  theory, 
especially,  in  generating  kilometer-scale  irregularities. 
The  measured  threshold  powers  for  exciting  a  few 
kilometer-scale  irregularities  in  the  Arecibo  experi¬ 
ments  agree  reasonably  well  with  the  theoretical 
values  of  Cragin  et  al.  However,  the  observed 
hundreds  of  meter-scale  irregularities  cannot  be  ex¬ 
cited  by  Cragin  et  al.  instability  because  the  required 
threshold  powers  Plh  are  too  large  according  to  the 
scaling  law:  Plh  x  /.L  *  (in  contrast  to  t.L  1  of  Perkins 
and  Valeo  [1974]).  The  detection  of  hundreds  meter- 
scale  irregularities  by  UHF  scintillation  technique 
was  recently  performed  by  Frey  et  al.  [1984]  in 
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Tromso  heating  experiments.  Frey  et  al.  considered 
three  candidate  mechanisms  including  Perkins  and 
Valeo  [1974],  Cragin  et  al.  [1977],  and  ours  [Kuo 
and  Lee ,  1983].  They  reported  that  the  thresholds  of 
our  thermal  filamentation  instability  were  consistent 
with  their  measured  ones,  whereas  those  of  the  other 
two  instabilities  were  higher  by  at  least  a  factor  of 
100 

It  has  been  predicted  that  the  conceptualized  Solar 
Power  Satellite  (SPS)  will  produce  large-scale  iono¬ 
spheric  irregularities  [Perkins  and  Goldman ,  1981] 
and  also  possibly  together  with  the  simultaneous  ex¬ 
citation  of  earth’s  magnetic  field  fluctuations  [ Lee 
and  Kuo.  1984/)].  A  rocket  experiment  named  Micro- 
wave  Ionosphere  Nonlinear  Interaction  Experiment 
(MINIX)  was  carried  out  recently  by  the  Kyoto  Uni¬ 
versity  group  in  Japan  to  simulate  the  SPS  [ Matsu - 
morn  et  al..  1984],  The  microwaves  were  transmitted 
at  a  frequency  of  2.45  GHz  with  incident  power  den¬ 
sities  comparable  to  the  envisioned  intensity  of  230 
W  m2  in  SPS.  The  beam  width  of  order  of  a  few 
meters  is,  however,  less  than  that  (  ~  10  km)  proposed 
for  SPS  by.  at  least,  3  orders  of  magnitude.  The  exci¬ 
tation  of  electron  cyclotron  waves  and  electron 
plasma  waves  were  observed,  but  it  is  not  surprising 
that  ionospheric  irregularities  and  other  predicted 
phenomena  could  not  be  produced  in  MINIX  be¬ 
cause  of  the  rather  narrow  microwave  beam.  The 
beam  size  is  probably  a  very  crucial  parameter  in  the 
simulation  experiments  of  SPS. 

4  CONCLUSIONS 

Thermal  filamentation  instability  can  be  excited  in 
space  plasmas  by  powerful  radio  waves  whose  fre¬ 
quencies  exceec.  half  the  electron  gyrofrequency. 
Large-scale  magnetic  and  plasma  density  fluctuations 
can  be  simultaneously  created  if  the  transmitter  of 
radio  waves  is  operated  continuously  for  a  few  min¬ 
utes.  During  its  linear  stage,  this  instability  is  able  to 
produce  significantly  large  earth's  magnetic  field  fluc¬ 
tuations  in  the  ionosphere  and  magnetosphere  that  is 
a  low -/(  plasma.  A  threshold  is  required  for  the  insta¬ 
bility  fundamentally  because  the  high-frequency  side¬ 
bands  and  the  purely  growing  modes  are  not  eigen- 
modes  rather  the  nonlinearly  driven  modes.  The  off- 
resonance  effect  (i.e.,  detuning)  imposes  the  primary 
damping  mechanism  of  the  instability.  The  electron- 
ion  (electron-neutral)  collisions  play  the  following 
roles  in  the  excitation  of  thermal  filamentation  insta¬ 
bility:  (1)  to  convert  the  radio  wave  energy  into  the 
driv  ing  heat  source  of  the  instability,  (2)  to  reduce  the 


electron  temperature  perturbations  induced  by  the 
heat  source:  this  is  also  the  principal  process  of 
energy  loss  in  determining  the  threshold  level  of 
large-scale,  highly  field  aligned  instability  (see  equa¬ 
tion  (3)),  and  it  is  therefore  commonly  termed 
“electron-ion  (electron-neutral)  cooling”  [e.g.,  Farley 
et  al..  1983],  and  (3)  to  facilitate  cross-field  mobility 
of  charged  particles  for  establishing  the  collective  os¬ 
cillations  (i.e.,  the  field-aligned  purely  growing 
modes)  and,  as  a  consequence,  the  growth  rate  of  the 
instability  is  generally  proportional  linearly  to  the 
electron-ion  (electron-neutral)  collision  frequency. 

We  have  proposed  the  thermal  filamentation  insta¬ 
bility  as  the  cause  of  HF  wave-induced  geomagnetic 
field  perturbation  in  Tromso  ionospheric  heating  ex¬ 
periments  [Stubbe  and  Kopka.  1981;  Stubbe  et  al., 
1982].  The  thresholds  of  this  instability  were  found 
to  agree  with  the  measured  ones  in  recent  Tromso 
experiments  [Frey  et  at..  1984],  However,  our  pro¬ 
posed  mechanism  and  the  one  suggested  by  Fejer  and 
Krenzien  [1982]  among  others  need  to  be  further 
verified  in  the  future  experiments  with,  for  instance, 
coordinated  in  situ  measurements  in  the  heated  iono¬ 
spheric  region.  The  thermal  filamentation  instability 
of  VLF  waves  (whistlers)  can  occur  in  the  mag¬ 
netosphere  with  growth  times  of  a  few  minutes.  Our 
work  thus  suggests  a  wave  plasma  interaction  pro¬ 
cess  for  future  VLF  wave  injection  experiments  to 
explore  with  the  operation  of  transmitters  in  a  CW 
rather  than  pulsed  wave  mode.  The  Japanese  simula¬ 
tion  experiments  of  Solar  Power  Satellite  [Matsu- 
moto  et  al..  1984]  unfortunately  cannot  test  the  Per¬ 
kins  and  Goldman's  predictions  and  ours  of 
microwave-induced  ionospheric  disturbances  because 
of  the  intrinsic  limitation  on  the  transmitter  size  in  a 
rocket.  The  space  station  under  consideration  by 
NASA  probably  can  provide  better  experimental 
conditions  for  the  future  simulation  experiments. 
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The  theory  of  resonant  electron  diffusion  as  an  effective  saturation  process  of  the 
auroral  kilometric  radiation  has  been  formulated.  The  auroral  kilometric 
radiation  is  assumed  to  be  amplified  by  the  synchrotron  maser  instability  that  is 
driven  by  an  electron  distribution  of  the  loss-cone  type.  The  calculated  intensity 
of  the  saturated  radiation  is  found  to  have  a  significantly  lower  value  in  com¬ 
parison  with  that  caused  by  the  quasi-linear  diffusion  process  as  an  alternative 
saturation  process.  This  indicates  that  resonant  electron  diffusion  dominates  over 
quasi-linear  diffusion  in  saturating  the  synchrotron  maser  instability. 


1.  Introduction 

Coherent  electromagnetic  radiation  can  be  generated  in  a  magnetized  non- 
Maxwellian  plasma  by  the  synchrotron  maser  instability.  This  instability  occurs 
under  the  conditions  of  inverted  population,  and  stems  from  the  relativistic  mass 
dependence  of  the  electron  gyrofrequency  that  gives  rise  to  phase  bunching  in  the 
electron  gyration  orbits  (see,  forexample,  Bekefi,  Hirshfield  &  Brown  1961).  The 
operation  of  this  instability  has  been  suggested  by,  for  instance,  Melrose  (1976) 
and  Wu  &  Lee  (1979)  to  be  responsible  for  auroral  kilometric  radiation  (see,  for 
example,  Gurnett  1974;  Benson  1985  and  references  therein). 

The  optimum  environment  for  the  excitation  of  this  instability  in  space 
plasmas  is  provided  by  the  following  processes  that  invert  the  population  and 
form  the  non-Maxwellian  electron  distribution  for  lasing.  Kilovolt  electric 
potential  drops  are  believed  to  exist  along  the  auroral  field  lines  during  the 
inverted-V  events  (e  g.  Gurnett  &  Frank  1973;  Lin  &  Hoffman  1 982  and  references 
therein).  Although  the  depletion  of  the  background  electrons  to  a  very  low  level 
is  expected  because  of  the  parallel  electric  field,  the  upgoing  electrons  that 
originate  in  the  reflected  plasma-sheet  electrons  introduce  a  loss-cone  distri¬ 
bution.  The  non-Maxwellian  electron  distribution  thus  created  is  capable  of 
producing  intense  coherent  electromagnetic  radiation  in  the  upper  atmosphere 
(Benson  1985). 


The  synchrotron  maser  instability  converts  the  electron  energy  into  coherent 
radiation.  More  specifically,  the  amplification  of  electromagnetic  radiation  occurs 
at  the  expense  of  the  kinetic  energy  of  electrons  whose  perpendicular  momentum 
falls  in  the  range  where  the  electron  distribution  function  has  a  positive  slope. 
The  saturation  of  this  instability  is,  therefore,  carried  out  by  processes  that  can 
decrease  the  positive  slope  of  the  momentum  distribution.  A  quasi -linear 
diffusion  process  has  been  analysed  by  Wu  et  al.  (1981)  toevaluate  the  saturation 
level  of  the  auroral  kilometric  radiation.  This  diffusion  process  suppresses  the 
instability  by  driving  electrons  into  the  lower  velocity  region,  i.e.  by  reducing 
the  positive  gradient  of  the  electron  distribution  function.  Since  the  radiation 
frequency  is  less  than,  but  very  close  to,  the  electron  gyrofrequency,  resonance 
broadening  was  also  examined  by  Wu  et  al.  as  another  potential  saturation 
process.  They  found  that  the  saturation  wave  energy  would  be  one  order  of 
magnitude  higher  if  resonance  broadening,  rather  than  quasi-linear  diffusion, 
were  considered  to  be  the  primary  saturation  process.  In  other  words,  quasi-linear 
diffusion  dominates  over  resonance  broadening  as  a  saturation  process  of  the 
synchrotron  maser  instability  in  generating  auroral  kilometric  radiation. 

A  computer  simulation  of  auroral  kilometric  radiation  was  performed  later 
by  Wagner  et  al.  ( 1 983).  It  shows  that  the  synchrotron  maser  instability  saturates 
when  the  resonant  electrons  diffuse  into  the  loss  cone  as  a  result  of  the  turbulent 
scattering  of  particles  by  the  amplified  radiation.  The  turbulent  scattering  means 
the  resonant  broadening  due  to  the  perturbations  on  the  particle  trajectories 
caused  by  the  amplified  radiation  whose  saturation  intensities  can  be  deter¬ 
mined  by  computer  simulation  and  analytically  from  our  equation  (17).  The 
resonant  electrons  referred  to  are  those  that  satisfy  the  synchrotron  maser 
resonance  condition:  cur~  Q.-k.vt  =  0  where  «r  is  the  real  part  of  the  radiation 
frequency,  Q  defined  by  eBJymc  is  the  relativistic  electron  gyrofrequency,  k 
and  v  represent  the  wavenumber  of  radiation  and  the  electron  velocity,  respec¬ 
tively.  The  subscript  z  denotes  the  z  axis  of  a  rectangular  reference  frame  taken 
along  the  geomagnetic  field.  The  parameters,  e,  m,  c,  B0,  and  y  defined  by 
(l-t>2/c*)-I,  have  their  respective  conventional  meanings  of  electron  charge, 
electron  rest  mass,  the  speed  of  light  in  vacuum,  the  geometric  field  intensity, 
and  the  relativistic  factor.  In  the  present  paper,  we  investigate  analytical!}'  this 
computer  simulated  physical  process.  An  upper  bound  of  the  corresponding 
saturated  radiation  can  then  be  derived  analytically,  and  compared  quanti¬ 
tatively  with  the  saturated  radiation  deduced  from  the  quasi-linear  diffusion 
saturation  process. 

The  formulation  of  the  resonant  electron  diffusion  is  presented  in  §2,  starting 
with  the  relativistic  Vlasov  equation.  This  is  followed  by  a  quantitative  analysis 
and  comparison  of  resonant  electron  diffusion  with  quasi-linear  diffusion  as  the 
dominant  process  that  causes  the  saturation  of  the  synchrotron  maser  instability 
in  the  production  of  auroral  kilometric  radiation.  A  summary  and  conclusions 
are  presented  in  §3. 
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2.  Formalism 

2.1.  Resonant  electron  diffusion 

Resonant  diffusion  of  the  background  electron  distribution  in  velocity  space 
results  from  the  interaction  of  electrons  with  the  amplified  electromagnetic 
radiation.  We  use  the  following  relativistic  Vlasov  equation 

|^  +  (p/ywi).  V-e[E  +  (p/ymc)  x  (B  +  zP0)].  VFj/  =  0,  (1) 

where  p  is  the  electron  momentum  and  B0  is  the  background  geomagnetic  field 
assumed  to  be  uniform.  E  and  B  designate  the  wave  fields  of  the  amplified 
radiation,  which  is  assumed  to  have  a  right-hand  circular  polarization,  that  is, 

E  *2(i  +  ip)ee*'*,-u>, 

*  (2) 

B  =  -  i  £  (A  +  i0)  (fcc/w)  e 

k 

where  u>{  =  u>(k))  and  e(  =  e(k))  are  the  wave  frequency  and  field  intensity  of  the 
amplified  radiation,  respectively. 

The  distribution  function  (/)  is  assumed  to  be  a  linear  combination  of  the 
background  ((f))  and  the  perturbation  (Sf)  portions: 

/=  </>+<v =/«+/i+<y0+<*A 

where  </)  =  /0  +fx  and  Sf  =  Sf0  +  Sfx.  The  four  parts  of  the  distribution  function, 
/„,  fv  Sf0  and  Sflt  correspond  to  the  unperturbed  background  distribution,  the 
perturbation  on  the  background  distribution  (a  second-order  effect),  the  phase 
coherent  linear  density  perturbation  (a  first-order  effect),  and  the  phase  coherent 
nonlinear  density  perturbation  (a  third-order  effect),  respectively.  The  un¬ 
perturbed  background  distribution  (/„)  is  of  the  loss-cone  type  given  by 

fa(PL,  P.)  =  ^n-i^P^^Pl1  exp  ( -  PVAPf),  (3) 

where  j(>  0)  is  the  loss-cone  index,  APj  =  [mTt/(\  +$j)]f,  and  Pz  = 

The  phase  coherent  linear  (£/„)  and  nonlinear  (<$/,)  density  perturbations  are 
caused  by  the  instability  via  wave-particle  resonant  interaction.  The  phase 
incoherent  density  perturbation  due  to  nonlinear  wave-wave  interaction  will 
not  be  included  in  our  analysis,  for  simplicity.  However,  it  should  be  pointed  out 
that  when  the  amplified  wave  field  grows  to  a  significantly  high  level,  the  wave- 
wave  interaction  process  may  also  become  important  and  contribute  additively 
to  the  saturation  of  the  instability. 

Along  the  unperturbed  trajectory  dr/dt  =  P /yrn  and  dP/dt  =  -  QP  x  z,  the 
governing  equations  deduced  from  (1)  for  perturbations  of  the  distribution 
function  are  given  by 


* 

(4) 

d  , 

dj '  ~  “ 

([£,  +  l 

(£c)H'v,<%+,/')+cc>' 

(5) 
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a*-‘[B+(£)xBW*’ 


(6) 


where  both  c.c.  and  the  superscript  *  indicate  the  complex  conjugate.  Integrating 
(4)  along  the  unperturbed  trajectory,  one  obtains  the  phase  coherent  density 
perturbation  due  to  the  instability, 

«*«*-*>  cie\L  kPA  £  kP±  <n 

*  k  (J  —  Q  —  kv1  ymoj)  BPX  ymoj8Ply° 

where  6  =  tan-1  { Py/Px ). 

Similarly,  integrating  (6)  leads  to 

[7  kPt\  8  kP±  d] 

1  ?  oj-£l-kvt  [\  ymoj)  3P±  ymojdp\1' 

Substituting  (7)  and  (8)  into  (6)  yields 

s'- "  f 


x  1  [1 

d 

i  _  _L 

.(JsL) 

8 

i 

(w-Cl-kv,) 

^  ymo>) 

8P± 

\ym<a} 

8Pt\ 

which,  after  being  integrated  along  the  unperturbed  trajectory,  becomes 

/, ,  _L+r/_i+.L\ 

h  f|»p  [p^pf^pyyWypJ 


</>•  (10) 


Substituting  (10)  into  (8),  one  finally  gets  the  phase  coherent  nonlinear  density 
perturbation  due  to  the  instability, 

Sf  -  y.  1*  -«„'_•»)<  el{k*~~t)  ig  T  j  8  j£M 

w-n-fcr, 

x  + [(sfe+?;) ("' + Q) + 

*  1111 


where  oj‘  ( =  ai'(k'))  and  e'  ( =  e'(k'))  are  the  wave  frequency  and  field  intensity  of 
radiation  with  wavenumber  k'.  One  can  note  from  (7)  and  (11)  that  the  wave- 
particle  resonant  interaction  occurs  at  wT  -  0  -  kv,  -  0  where  air  is  the  real  part 
of  the  wave  frequency  oj. 


2.2.  Saturation  of  the  instability 

The  linear  and  nonlinear  density  perturbations,  <J/0  and  Sfv  give  rise  to  the 
induced  current  density  defined  by 


SJ  = 


dOPJP/ym)6f, 
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where  »0  is  the  unperturbed  background  electron  density.  We  shall  show  that 
the  nonlinearly  induced  current,  imposing  a  stabilizing  effect,  determines  the 
saturation  level  of  the  instability.  These  induced  currents,  appearing  in  the 
following  wave  equation,  acts  as  the  driving  source  of  the  radiation: 

jR-enro  - (12) 


With  the  substitution  of  the  calculated  induced  current  densities,  this  wave 
equation  yields  the  nonlinear  dispersion  relation  for  the  synchrotron  maser 
instability.  It  is  found  to  be 

-■-ft. -  ^4hdP-y 

+ ep, k ”■]  K  -  a  -  ft,|’  [Q  iPi + w ■  IP,]!  <^>)  ]  ,13) 

where  w  =  w,  +»r,  the  real  and  imaginary  parts  being  the  real  wave  frequency 
and  the  nonlinear  growth  rate,  respectively. 

Solving  (13)  for  wr  and  f,  we  derive  the  nonlinear  growth  rate: 


P  _  " 


S(P±-P±.) 


where 


+ [(ir- +  r)  1 2-;-  2f».  -  at + gs  v.,] 

'  (adPL  +  k  • 


8P , 


(14) 


lcc  +  wr(  Qjj  +  kV  -  u>,))  i 


-«c[ 

[  \wr  mu)rJ  mV] 


]• 


The  following  identities  have  been  used  in  the  derivation  of  (14): 
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and 


for 

where  P  means  the  principal  part.  We  remark  that  the  real  wave  frequency  a>r 
of  the  amplified  radiation  is  determined  primarily  by  the  bulk  electrons  of  the 
distribution  rather  than  the  resonant  electrons.  It  is  then  reasonable  to  assume 
that  the  real  wave  frequency  is  only  weakly  affected  by  the  nonlinearity  due 
to  the  wave-particle  interaction.  Therefore,  we  will  not  include  the  effect  of 
nonlinear  frequency  shift  in  the  present  work  and  simply  approximate  cur  with 
the  value  derived  from  the  linear  dispersion  relation  deduced  from  (13)  after 
dropping  the  nonlinear  terms. 

The  saturation  condition  for  the  instability  is  determined  by  (14)  after  setting 
T  =  0.  The  result  has  an  unknown  integral  over  the  average  background  distri¬ 
bution  function,  (/).  An  estimate  of  it  can  be  reasonably  made  on  the  basis  of 
the  following  arguments.  The  growth  rate  of  the  instability  and  the  distribution 
function  can  be  assumed  to  re  main  unchanged  for  most  of  the  instability  evolution 
time.  It  is  only  during  the  final  stage  of  the  instability  that  the  growth  rate  drops 
sharply  to  zero  with  the  formation  of  a  plateau  in  the  velocity  range  of  the 
distribution  function  where  a  positive  slope  existed  originally.  Mathematically 
while  | a)  —  Q  -  kvz |2  ->■  0  because  of  the  resonance  condition  and  P  ->•  0  at  the  final 
stage  of  the  instability,  [iYd/cP x  +  kvx8/8Pz]  (f)  -*■  0  because  of  the  formation  of 
a  plateau  in  the  distribution  function.  For  an  order  of  magnitude  estimate  of 
the  saturation  level,  we  let 


xi..  r>  ; _ \  _  y[y  Qq/ a)T  + kPz/mu)^  ^  D  n 

-  Q - kvz)  - - - - S{P±  - P±0) 


p,>  <  p,  <  p„. 


K-Q- k'vt\* (D dPx  +  k  *’x dp)  -  r?  (Q 8PX  +  k  v± ep)f°  ( 1 5) 


in  (14),  where 


fo  (16) 


is  the  initial  (namely,  the  linear)  growth  rate  of  the  instability.  The  assumption 
that  the  growth  rate  of  the  synchrotron  maser  instability  driven  by  a  loss-cone 
distribution  remains  at  its  initial  value  for  a  large  fraction  of  the  time  has  been 
justified  in  Wu  et  al.  (1981). 

After  a  few  exponentiation  periods,  the  fastest  growing  wave  becomes 
dominant.  Since  the  frequency  spectrum  of  the  amplified  radiation  narrows  down 
sufficiently,  the  single-wave  approximation  is  good  enough  for  the  purpose  of 
estimating  the  saturation  level  of  the  instability.  In  view  of  those  points  and 
noting  that  T ,(k)  is  largest  at  k  =  0  (Wu  1981),  we  can  further  simplify  the 
analysis  of  (14)  by  considering  only  a  single  wave  with  1=0  for  estimating  the 
upper  bound  of  the  saturated  wave  intensity. 
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2.3.  Intensity  of  the  saturated  radiation 

With  these  simplifications  the  integrations  in  (14)  can  be  carried  out  analytically. 
Then,  taking  T  =  0  in  (14)  leads  to  the  following  expression  for  the  saturation 
intensity: 


where 

and 


e! 


4irn0Te 


5*3[j-2««(j+l)/(2j  +  3)]|  5j?r  2(j+l)(3j  +  2)  , 
<o%  2j(2j  +  3)(3j+l)  \  y?l  j(2j+  3)  (3j  +  1)“ 


4(j+l)U  +  2)a*  (j—l)(2j+ 1) 

j(2j  +  3)(2j+5)(3j+l)  2(3 \j  +  1 )  a* 


O'*1)*2  (: 

yf  \  <%Jj(2j+3)(3j+t)\J 


—  a* 


s±4\p 

2j  +  5/1 


*V  =  wr/Q0,  Wp  =  Wp/Q0)  y,  =  rya0)  /?2  =  kP^/mh* 


a*=(l-S*)//?*<5*. 


(17) 


The  normalized  linear  growth  rate  y,,  obtained  from  (16),  is  given  by 


yi  =  L-_a.^±_2l 

r'  /?*wr  (2j  + 1) !  l  2jf  +  3  J 


(18) 


It  is  clear  from  (18)  that  y,  is  negative  for  j  =  0  corresponding  to  the  case  of  a 
Maxwellian  distribution  (see  (2)).  The  positive  index  j  characterizes  the  velocity 
spread  of  the  loss-cone  distribution  (defined  by  (2))  in  the  v±  space.  It  is  the 
population  inversion  in  Px  that  provides  the  source  of  free  energy  for  the  syn¬ 
chrotron  maser  instability. 

We  now  evaluate  the  intensity  of  the  saturated  radiation  based  on  (17)  with 
the  intention  of  comparing  resonant  diffusion  with  quasi-Iinear  diffusion  as  the 
dominant  saturation  process  in  auroral  kilometric  radiation  that  is  amplified  by 
the  synchrotron  maser  instability.  The  typical  magnetospheric  parameters  that 
were  used  by  Wu  (1981)  to  analyse  quasi-Iinear  diffusion  in  the  loss-cone  distri¬ 
bution  with  j  =  1  are  also  adopted  as  follows: 

<op  —  0-2,  Tt  =  5keV,  n0  =  1  cm-*,  ft1  ~  0  02. 

For  the  k  =  0  mode,  the  calculated  3r  is  about  0-999  and  then  a  ~  0-112  and 
y(  ~  6-8x  10~3.  With  these  parameters,  we  obtain  the  following  results  from  (17): 

.1  J 

7 — '—=  -  8-7  x  10-6  or  ~  7-0x  10-13(ergcm~3)  (19) 

4nnaT.  4  n  ° 


for  n0  =  1  cm-3  and  Tt  =  6  keV. 

This  energy  density  (7-0  x  10~l3ergcm-3)  of  the  saturated  radiation  is  lower 
than  those  ( ~  10-10ergcm_3)  obtained  from  the  consideration  of  the  quasi- 
Iinear  diffusion  process  by  about  two  orders  of  magnitude.  This  indicates  that 
resonant  diffusion  is  more  effective  than  quasi-Iinear  diffusion  in  causing  the 
saturation  of  the  synchrotron  maser  instability  that  is  driven  by  a  loss-cone 
distribution. 
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3.  Summary  and  Conclusions 

We  have  formulated  the  theory  of  resonant  electron  diffusion  as  a  saturation 
process  of  the  synchrotron  maser  instability  driven  by  a  loss-cone  electron  distri¬ 
bution.  We  were  motivated  by  the  computer  simulation  work  of  Wagner  el  al. 
(1983).  They  showed  that  auroral  kilometric  radiation  amplified  by  the  syn¬ 
chrotron  maser  instability  saturates  when  the  resonant  electrons  diffuse  into  the 
loss  cone  via  turbulent  scattering  of  electrons  by  the  excited  radiation.  We  are 
able  to  evaluate  analytically  the  intensity  of  the  saturated  auroral  kilometric 
radiation.  The  ratio  (~  8-7  x  10-s)  of  the  saturated  wave  energy  to  the  kinetic 
energy,  as  given  in  (19),  is  less  than  that  (~  20x  10-3)  obtained  in  computer 
simulation  by  a  factor  of  about  20  because  different  distribution  functions  and 
parameters  are  used. 

The  main  purpose  of  this  paper  is  to  find  out  whether  resonant  diffusion  indeed 
dominates  over  quasi-linear  diffusion  in  causing  the  saturation  of  auroral  kilo¬ 
metric  radiation.  The  saturated  radiation  intensity  has  been  calculated  with  the 
same  distribution  function  and  plasma  parameters  that  wrere  employed  in  Wu 
(1981)  for  the  analysis  of  the  quasi-linear  diffusion  process.  The  radiation 
saturated  by  the  resonant  diffusion  process  has  a  significantly  lower  intensity 
in  comparison  with  that  by  the  quasi-linear  diffusion  process.  We  conclude, 
therefore,  that  resonant  electron  diffusion  provides  the  dominant  process  for 
saturating  the  synchrotron  maser  instability  in  the  generation  of  auroral  kilo¬ 
metric  radiation  as  Wagner  et  al.  (1983)  demonstrated  in  the  computer  simu¬ 
lation. 
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